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Summary: 

This  project  has  resulted  in  a  solid  improvement  to  the  state  of  the  art  knowledge  in 
nonlinear  structural  dynamics.  The  new  insight  gained,  and  the  types  of  structural 
components  analyzed,  are  directly  applicable  to  the  structural  mechanics  objectives  of  the 
Air  Force.  In  order  to  fully  provide  a  test-bed  for  future  hypersonics  vehicular  analysis 
(in  particular)  it  is  crucial  to  understand  the  nonlinear  dynamical  behavior  of  slender 
structural  elements  in  extreme  environments. 

The  technical  papers  that  form  the  appendices  of  this  final  report  detail  the  achievements 
that  have  been  made  on  both  the  experimental  and  computational  fronts  of  this  project. 
Both  Pi’s  have  spent  extended  periods  at  the  AFRL  in  Dayton,  Ohio  interacting  with  Air 
Force  personnel  in  this  area  (including  during  the  summers  of  2009,  2010  and  2011).  A 
number  of  research  assistants  were  supported  on  this  project.  Both  Pi’s  and  the  two  of  the 
graduate  students  spent  time  at  AFRL  in  Dayton  in  the. 

The  Context  of  the  Research: 

In  order  to  meet  ever-increasing  performance  demands  modern  aircraft  are  being 
designed  lighter  than  ever  before.  This  weight  efficiency  usually  comes  at  the  expense  of 
structural  mass,  leading  to  more  slender  structural  components.  Slender  structures  are 
however  susceptible  to  vibration  and  instabilities,  particularly  buckling,  which  can  occur 
well  before  strength  limits  are  reached. 

Snap-through  buckling  is  a  particular  type  of  buckling  where  a  structure  snaps  from  one 
state  to  another  remote  state.  Such  large  deflections  pose  a  hazard  since  they  can  cause 
fatigue,  and  in  the  context  of  hypersonics,  sonic  fatigue  due  to  boundary  layer  effects. 
Postbuckled  snap-through  occurs  in  axially  loaded  structures,  where  a  structure  snaps 
between  the  two  (often)  symmetric  buckled  states  under  some  perturbation. 

The  two  main  thrusts  of  this  research  has  concerned  low-order  experimental  models,  and 
high-order  computational.  Although  the  approaches  and  scale  of  these  system  approaches 


are  quite  different,  the  resulting  research  discoveries  illustrate  the  similarity  in  dynamic 
phenomena,  and  indeed,  these  approaches  have  mutually  informed  each  other. 

Summary  of  Experimental  Results: 

Figure  1(a)  shows  a  photograph  image  of  the  experimental  test-bed,  built  by  the  PI  Prof. 
Virgin  and  his  students  at  Duke  University.  A  schematic  of  this  system  is  shown  in  part 
(b).  When  this  system  is  laterally  excited,  it  is  capable  of  exhibiting  a  sensitivity  to  initial 
conditions.  Part  (c)  shows  how  the  system  can  be  randomly  perturbed  about,  and 
between,  various  co-existing  behaviors.  That  is,  the  system  ‘snaps’  between  various 
forms  of  behavior,  and  this  of  course  is  exactly  the  type  of  behavior  that  has  such  striking 
implications  for  sonic  fatigue. 


Figure  1:  The  link  model,  (a)  the  system  as  constructed  in  the  lab,  (b)  a  schematic,  (c)  a 
time  series  in  which  the  system  is  perturbed  between  co-existing  states. 

More  details  of  the  types  of  dynamic  behavior  relating  to  this  system  can  be  found  in 
references  [1-4]. 


Summary  of  Numerical  Results 


A  sophisticated  finite  element  analysis  and  simulations  of  snap-through  behavior  was 
conducted  by  Prof.  Stanciulescu  (the  subcontractor)  and  her  students  at  Rice  University. 
In  addition  to  the  mechanical  instability  of  interest  progress  was  also  made  on  accounting 
for  the  coupling  between  mechanical  and  thermal  loading. 

As  an  example,  Figure  2(a)  shows  through  an  indirect  evaluation,  a  lower  bound  on  the 
(snap-through)  temperature  that  can  be  identified  for  a  given  shallow  arch.  This  is  the 
continuous  analog  of  the  discrete  system  discussed  in  the  experimental  section.  Below 
these  temperatures  the  arch  will  not  experience  snap-through.  A  direct  method  for  tracing 
the  stability  boundary  of  coupled  systems  was  also  formulated  and  implemented.  It  can 
handle  any  type  of  critical  points  and  also  identifies  mode  changes,  with  a  typical  set  of 
results  displayed  in  Figure  2(b). 


Figure  2:  Sample  finite  element  results,  (a)  Safe  snap-through  boundaries  for  a  shallow 
arch  including  thermal  loading,  (c)  Critical  loads  under  various  conditions. 

Summary  of  Experimental/Numerical  Correlation 

One  of  the  original  goals  of  this  research  project  was  to  utilize  parallel  approaches  to 
experimental  and  numerical  studies.  A  typical  example  is  shown  in  Figure  3.  The  lower 
data  (plotted  in  black)  shows  how  the  response  of  the  experimental  system  depends 
strongly  on  the  forcing  parameters.  The  response  of  the  system  is  monitored  in  terms  of 
the  number  of  ‘snaps’  per  forcing  cycle.  For  example,  when  the  system  is  forced  at  a 
fixed  amplitude,  at  a  frequency  £2  of  4  Hz,  no  snapping  behavior  occurs  and  the  system 
simply  oscillates  (with  small  amplitude)  about  one  of  the  (buckled)  equilibria.  Over  the 
forcing  frequency  range  £2  =  4.4  to  6.2  the  system  may  either  not  snap-through,  or  snap- 
through  periodically  (hence  the  ‘0’  or  ‘2’  label).  However,  over  the  range  £2  =  6.3  to  9 
Hz, 


The  system  exhibit  very  complex  (sometimes  chaotic)  behavior.  Some  sample  time  series 
are  shown  as  insets  in  this  figure.  The  upper  part  of  the  diagram  shows  the  corresponding 
numerical  result.  The  correlation  is  remarkable. 
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Figure  3:  A  comparison  between  experimental  and  numerical  results. 

Summary: 

The  research  achieved  in  this  program  has  resulted  in  a  range  of  innovative  approaches  in 
both  the  experimental  and  computational  realms.  This  work  provides  a  strong  foundation 
on  which  to  build  more  sophisticated  studies  based  on  more  realistic  (panel  and  shell-like 
structures,  including  composites).  The  list  of  publications  from  this  research  are  given 
below,  and  be  consulted  for  further  details. 

Publication  List: 


(1)  R.  Wiebe,  L.N.  Virgin,  I.  Stanciulescu,  and  S.M.  Spottswood,  ‘On  Snap-Through 
Buckling’,  52nd  AIAA/ASME/ASCE/AHS/ASC  Structures,  Structural  Dynamics 
and  Materials  Conference,  Denver,  CO  (2011). 

(2)  Y.  Chandra,  I.  Stanciulescu,  T.G.  Eason,  and  S.M.  Spottswood,  ‘Numerical 
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504,2011. 
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loads’ ,  International  Journal  for  Numerical  Methods  in  Engineering ,  accepted  for 
publication. 
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ABSTRACT 

Geometrically  nonlinear  structures  often  possess  multiple  equilibrium  configurations .  Under  extreme  condi¬ 
tions  of  excitation  it  is  possible  for  these  structures  to  exhibit  oscillations  about  and  between  these  co-existing 
configurations.  This  behavior  may  have  serious  implications  for  fatigue  in  the  context  of  aircraft  surface  panels. 
Snap-through  is  a  name  often  given  to  sudden  changes  in  dynamic  behavior  associated  with  mechanical  instability 
(buckling).  This  is  an  often  encountered  problem  in  hypersonic  vehicles  in  which  severe  thermal  loading  and  acous¬ 
tic  excitation  conspire  to  create  an  especially  hostile  environment  for  structural  elements.  In  this  paper ,  a  simple 
link  model  is  used ,  experimentally  and  numerically,  to  investigate  the  mechanisms  of  snap-through  buckling  from  a 
phenomenological  standpoint . 


1  Introduction 


In  order  to  meet  ever  increasing  performance  demands  modern  aircraft  are  being  designed  lighter  than  ever  before.  This 
weight  efficiency  usually  leads  to  more  slender  structural  components.  Slender  structures  are  susceptible  to  instabilities, 
particularly  buckling,  which  can  occur  well  before  strength  limits  are  reached.  The  nonlinear  vibration  of  buckled  structures 
may  also  lead  to  a  large  variety  of  chaotic  and  periodic  oscillations.  We  are  primarily  interested  in  the  circumstances  under 
which  a  curved  structure  is  forced  to  its  inverted  conhguration,  as  shown  in  Fig.  1(a)  and  (b). 

Snap-through  buckling  is  a  particular  type  of  instability  where  the  structure  snaps  from  one  state  to  another  (remote) 
state.  Such  large  deflections  pose  a  hazard  since  they  can  exacerbate  fatigue  [1,2].  The  two  primary  types  of  snap-through 
buckling,  post  buckled/bifurcated  snap-through  and  limit  point  snap-through,  are  shown  in  parts  (c)  and  (d)  of  Fig.  1  [3] . 

Limit  point  buckling  occurs  in  structures  where  the  stiffness  decreases  (with  increasing  loading)  to  a  vanishing  point, 
i.e.,  a  horizontal  tangency  in  the  force-deflection  relation  (as  shown  in  Fig.  1(c)).  At  the  limit  point  the  structure  jumps 
to  the  remaining  stable  equilibrium.  At  loads  below  the  limit  point  the  structure  may  also  snap-through  under  external 
perturbations.  Postbuckled  snap-through  occurs  in  axially  loaded  structures,  where  a  structure  snaps  between  the  two  (often) 
symmetric  buckled  states  under  some  perturbation,  for  example,  from  point  be  to  6*  in  Fig.  1(d).  Postbuckled  snap-through 
typically  occurs  after  pitchfork  buckling;  and  it  often  requires  that  the  structure  first  be  loaded  to  a  buckled  conhguration  (as 
shown  in  Fig.  1(d)).  Related  behavior  can  also  occur  in  aircraft  surface  panels  that  may  buckle  due  to  thermal  loading,  [4-6] . 

When  a  structure  with  an  underlying  snap-through  static  behavior  is  subjected  to  excitation,  then  highly  nonlinear  oscil¬ 
lations  are  possible.  In  order  to  gain  phenomenological  insight  into  this  type  of  response  and  enable  relatively  unambiguous 
experimental  verification  we  introduce  a  very  simple  single  degree  of  freedom  (SDOF)  discrete  structural  system  (Fig.  2). 
This  system,  for  a  certain  induced  pre-stress  through  the  axial  spring  2 ks,  and  a(t)  =  constant  (i.e.,  the  actuator  is  held  in  a 
fixed  position),  possesses  two  symmetric  static  equilibria,  0  =  ±0e.  Furthermore,  when  the  system  is  subject  to  ’lateral’  exci¬ 
tation  a(t),  we  observe  the  possibility  of  complicated  behavior  involving  motion  that  is  heavily  influenced  by  the  locations  of 
the  underlying  equilibria.  A  key  aspect  of  this  form  of  dynamic  buckling  is  the  existence  of  an  unstable  equilibrium  between 
the  stable  states.  These  unstable  equilibria  define  the  basins  of  attraction  for  the  stable  states.  Hence,  the  energies  of  the  stable 
and  unstable  states  may  be  useful  for  distinguishing  between  trajectories  that  snap  from  those  that  do  not.  The  static  equilib¬ 
ria,  both  stable  and  unstable,  are  therefore  of  great  importance  in  studying  the  systems’  forced  dynamics.  This  system  will 
be  scrutinized  using  classical  analysis  techniques  as  well  as  experimental  verification  to  shed  light  on  the  issue  of  dynamic 
snap-through  in  the  larger  context  of  aircraft  structural  components.  Later  work  will  focus  on  the  multi-degree-of-freedom 
system,  in  which  mode  shapes,  asymmetric  behavior,  etc.  will  be  important  features. 

Single  degree  of  freedom  systems,  particularly  ones  with  Duffing-type  potential  wells,  are  well  represented  in  the  lit¬ 
erature.  An  excellent  review  of  the  behavior  of  the  forced  Duffing  oscillator  is  available  in  [7,8],  The  route  to  potential 
well  escape  (i.e.,  snap-through  for  a  mechanical  system)  is  studied  in  [9]  using  Melnikov  theory.  Melnikov  theory  is  also 
used  in  [10]  to  give  a  criterion  for  the  existence  of  chaos  in  quasi-periodically  forced  Duffing  oscillator.  An  SDOF  structure 
with  an  underlying  Duffing  potential  well  was  studied  in  [11],  This  paper  discusses  the  use  of  homoclinic  orbits  as  the 
bifurcation  indicator  in  an  unforced  system,  and  bifurcation  to  and  from  chaos  with  respect  to  harmonic  forcing  parameters. 


A  (modified)  Duffing  oscillator  was  studied  numerically  in  [12, 13].  In  the  vicinity  of  basin  boundaries,  small  changes  in 
initial  conditions,  parameters,  or  indeed  in  the  numerical  algorithm  used,  may  lead  to  different  attractors.  A  study  of  the 
bifurcations  in  forcing  parameter  space,  similar  to  the  goal  of  this  research,  is  done  in  [14] .  This  work  produced  approximate 
analytical  boundaries  separating  bounded  from  unbounded  motion. 

2  Experimental  Setup 

Shallow  arches  typically  snap-through  in  a  symmetric  mode  shape  [15-19],  and  thus  a  SDOF  link  arch  may  be  used  to 
qualitatively  model  this  behavior.  For  deeper  arches,  i.e.,  arches  with  a  higher  rise-to-span  ratio,  an  asymmetric  mode  of 
buckling  is  typically  encountered.  Figure  3  is  a  photograph  of  the  experimental  SDOF  link  model  tested.  The  two  link  arms 
were  both  carbon  fiber  and  had  a  length  of  32.4  cm  and  a  mass  of  41g.  The  center  joint  had  a  mass  of  166g  and  the  slider 
had  a  mass  of  1260g.  Note  that  the  system  was  set  up  horizontally  on  a  test  bed  and  hence  force  due  to  gravity  was  out  of 
the  plane  of  the  system  dynamics. 

The  two  axial  springs  are  an  essential  component  of  the  structure  and  result  in  the  underlying  multivalued  potential 
energy  surface  which  mimics  a  shallow  arch.  As  a  practical  consideration  it  was  decided  to  use  two  opposing  mutually  pre¬ 
tensioned  springs  (each  with  ks  =  74  N/m)  to  avoid  inducing  compressive  forces,  which  would  of  course  buckle  the  springs. 
The  two  transverse  springs  connected  to  the  center  joint  provided  a  transmitted  force  to  the  system  via  the  displacement 
of  the  opposite  end  of  the  active  spring,  ka.  The  displacement  was  provided  by  a  Scotch-yoke  mechanism  moving  with  an 
approximately  sinusoidal  motion  [8],  The  second  (anchor)  spring,  kp ,  ensured  a  pre-tension  in  the  active  spring.  The  active 
and  anchor  springs  had  stiffnesses  of  ka  =  38  N/m  and  kp  =  23  N/m  respectively.  The  control  parameters  for  the  system  were 
the  Scotch-yoke  amplitude  and  frequency  which  equate  to  the  forcing  amplitude  and  forcing  frequency  of  the  system. 

The  most  natural  generalized  coordinate  for  this  system  is  the  angle  of  the  link  arm  on  the  non-sliding  side  of  the 
structure.  This  angle,  0,  was  measured  using  a  high  speed  camera  (Prosilica  GC640).  The  arm  position,  and  hence  the  angle, 
was  obtained  by  locating  a  small  white  target  on  the  link  arm  (see  Fig.  3)  by  its  contrast  with  the  background  using  Lab  VIEW 
software.  This  process  was  done  in  real-time  for  each  frame  taken  by  the  camera,  providing  a  time  series  of  the  arm  angle. 
Experimental  results  are  discussed  later  in  sections  4  and  5. 

3  Modeling 

A  numerical  model  was  developed  to  compare  with  the  experimental  results.  Referring  back  to  Fig.  2,  we  model  the 
discrete  dynamic  system  using  an  energy  approach.  Here,  mi,  is  the  mass  of  the  link  arms,  mc  is  the  mass  of  the  center  joint 
(assumed  to  be  a  point  mass),  ms  is  the  mass  of  the  slider  assembly  and  joint,  L  is  link  length,  ks  are  the  lateral  structural 
spring  stiffnesses,  ka  is  the  active  spring  stiffness,  kp  is  the  anchor  spring  stiffness,  and  0  is  the  angle  of  the  link  arm  on 
the  non-sliding  side.  The  spring  masses  were  considered  negligible.  The  displacement  at  the  end  of  the  active  spring  a(t) 
is  transmitted  as  a  force  through  the  active  spring  to  provide  the  system  forcing.  The  angle  0o  is  the  unforced  equilibrium 
angle,  i.e.,  the  angle  of  the  stable  equilibrium  with  the  forcing  in  the  neutral  position.  Once  again  the  effect  of  gravity  is 
ignored  in  the  analysis  since  the  experimental  system  was  laying  flat.  The  kinetic  and  potential  energies,  along  with  the 


energy  dissipation  of  the  system  are  given  by 


T(Q,Q,t)  =  jj  (mc  +  \rrib) T2  +  \msx!s (0)2 

+\mb  4(0)2+4(0K(0)+^(0)2+x'(0)2  }0x 

V (0, f; 0i )  =  ks[xs(Q)  -Xs(0i)]2  +  \ka[a{t)  -yc(0)]2  +  \kpyc{Q)2, 

F(0,0,O  =  |p[x,(0)]2  =  ip[^(0)0]2, 


(1) 


where  xc(0)  =  Lcos(0),  yc(0)  =  Lsin(0),  xs(0)  =  2Lcos(0).  The  masses  of  the  springs  were  small  and  their  contribution  to 
the  kinetic  energy  was  neglected .  The  slider  mechanism  was  by  far  the  largest  contributor  to  the  damping  of  the  experimental 
system;  therefore  the  dissipation  function  F  was  set  to  be  a  function  of  xs  only.  The  angle  0i  is  the  angle  at  which  the 
two  lateral  structural  springs  were  in  equilibrium.  This  is  the  angle  at  which  the  potential  energy  of  these  two  springs  is 
at  a  minimum,  i.e.,  the  stable  equilibrium  of  the  system  if  the  transverse  (forcing)  springs  were  removed.  However  the 
equilibrium  states  observed  experimentally  are  at  the  stationary  points  of  the  total  potential  energy  [20] .  The  two  angles  0o 
and  0i  are  related  by  the  expression 


3 
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[v  (0;0i)]0=0o  =  °- 


(2) 


Note  that  the  two  forcing  springs,  ka  and  kp,  were  set  to  be  in  equilibrium  for  0  =  0. 

The  potential  and  kinetic  energy  equations  may  be  substituted  into  the  Euler-Lagrange  equation  [21] 


d  (d(T-V)\  dT  dV  dF 
dt  V  /  30  +  30  + 


(3) 


to  yield  the  equation  of  motion  of  the  system: 


\mc  +  | mb  +  2 (mb  +  2 ms)  sin2  0]  0  +  {mi,  +  2 ms)  sin  (20)02  +  4p0  sin2  0 
+8&s(cos0i  — cos0)sin0+  \  {ka+kp)  sin20—  {ka/L)cosB)a(t)  =0. 


The  damping  coefficient  p  (Kg/s)  was  determined  by  fitting  a  simulated  (using  fourth-order  Runge-Kutta  time  stepping 
of  Eq.  (4)  with  At  =  0.001  s)  free  decay  response  with  the  experimental  response.  In  order  to  check  the  numerical  stability  of 
the  Runge-Kutta  routine  it  was  also  compared  with  simulations  done  using  the  NDSolve  function  available  in  Mathematica, 


with  excellent  agreement.  Figure  4  (a)  shows  this  comparison  between  simulated  and  experimental  time  series  (for  0o  =  26.0° 
and  0i  =  36.2°).  We  note  that  the  agreement  is  good,  despite  the  large-amplitude,  snap-through  characteristic  of  the  free 
decay.  Parts  (b)  and  (c)  of  this  figure  show  the  sensitivity  of  the  simulation  error  norm  (calculated  as  the  average  of  the 
absolute  difference  between  simulated  and  experimental  time  series  up  to  t  =  6  s)  as  a  function  of  (3.  Part  (b)  is  for  the 
large  amplitude  experimental  data  points  shown  in  part  (a).  The  large  jumps  in  the  error  occur  because  as  the  damping  is 
changed  the  system  will  snap-through  either  too  often  (damping  too  low),  or  not  enough  (damping  too  high).  The  lower 
error  region  near  |3  =  0.65  occurs  because  the  system  in  fact  snaps-through  an  extra  two  times  and  therefore  returns  to  the 
same  side  as  the  experimental  data  so  the  error  is  actually  lower  than  for  higher  values  of  |3  where  only  one  extra  snap 
event  occurs.  This  leads  to  a  smaller  error  norm  even  though  it  is  less  accurate.  Note  that  although  the  system  seems  to  be 
very  sensitive  to  the  damping  coefficient,  this  plot  corresponds  to  a  single  experimental  time  series.  Finally  part  (c)  shows 
a  similar  sensitivity  study  done  using  a  small  non-snapping  free  decay,  which  shows  a  clear  minimum.  Parts  (a)  and  (b) 
hint  that  an  optimal  damping  is  approximately  (3  =  1.2,  whereas  part  (c),  for  the  small  amplitude  decaying  oscillation,  has 
a  optimal  value  closer  to  |3  =  1.3.  This  indicates  that  the  system  damping  is  slightly  larger  for  small  oscillations.  However, 
since  the  results  to  be  shown  tended  to  be  dominated  by  large  amplitude  oscillations  the  value  (3  =  1.2  was  used  throughout. 
The  model  could  be  refined  by  including  the  effects  of  Coulomb  damping  such  as  in  [22] ,  where  Coulomb  friction  is  shown 
to  induce  follower  forces  on  the  system.  However,  due  to  the  overall  excellent  agreement  between  the  experimental  and 
numerical  results  shown  later,  it  was  decided  that  the  viscous  friction  used  in  the  model  is  accurately  portraying  the  system 
dynamics  in  general,  including  energy  dissipation. 

4  Equilibria  and  Stability 

The  static  force-displacement  relationship  may  be  obtained  by  finding  the  static  force  (or  in  this  case  static  displacement 
a(t)  — >  ae)  and  angle  (),.  for  which  the  potential  energy  of  the  system  exhibits  a  stationary  point.  The  force-displacement 
relationship  is  given  by 
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(5) 


According  to  the  theorem  of  minimum  potential  energy  the  stability  of  the  equilibria  are  given  by  the  sign  of  the  second 
derivative  of  the  potential  energy  function  [20],  A  positive  (negative)  second  derivative  implies  a  minimum  (maximum) 
potential  and  hence  stable  (unstable)  equilibrium.  An  equivalent  approach  is  to  investigate  the  local  natural  frequency  of  the 
system  for  small  oscillations  about  the  equilibria  under  the  static  load  being  considered  [3,23],  This  method  is  preferred 
because  unlike  potential  energy  it  is  easy  to  directly  measure  the  response  frequency  of  an  experimental  system.  The  natural 
frequency  may  also  be  of  interest  in  its  own  right.  The  sign  of  to2  determines  stability,  i.e.,  if  to2  >  0  then  disturbances  do 
not  grow  in  time.  Linearizing  the  equation  of  motion  about  a  static  equilibrium  position  0  =  0^  under  a  static  load  a(t)  =  ae 
produces  the  equation  of  a  simple  harmonic  oscillator  from  which  the  undamped  natural  frequency  is  easily  determined  to 
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where  the  numerator  is  the  second  derivative  of  the  potential  energy.  Given  that  the  denominator  of  equation  (6),  i.e.,  the 
inertia,  is  positive  definite,  this  result  is  clearly  in  agreement  with  the  theorem  of  minimum  potential  energy. 

Figure  5  shows  the  relationship  between  the  force  (static  displacement  of  the  actuator  ae),  displacement,  and  local  natural 
frequency  for  both  the  experimental  system  and  the  analytical  model  (for  0o  =  26.0°  and  0i  =  36.2°).  The  experimental 
response  frequency  is  necessarily  the  damped  natural  frequency,  whereas  in  the  analytical  result  is  the  undamped  natural 
frequency.  However,  because  the  damping  is  relatively  small  these  two  values  are  still  very  close.  For  the  experimental 
results  there  are  no  data  points  in  the  unstable  region  given  that  these  equilibria  cannot  be  observed.  However  the  agreement 
is  excellent  in  the  stable  region. 

The  static  equilibria  of  dynamical  systems,  and  underlying  potential  energy,  give  much  insight  into  the  dynamic  re¬ 
sponse,  i.e.,  they  act  as  an  organizing  framework  for  the  nonlinear  behavior.  If  the  system  is  subject  to  periodic  excitation, 
they  also  influence  large  amplitude  snap-through  oscillations  and  highly  nonlinear  global  behavior.  This  will  be  considered 
in  the  next  section. 


5  Dynamic  Response 

The  forced  dynamic  response  of  the  system  was  simulated  by  applying  a  fourth-order  Runge-Kutta  time  stepping  scheme 
on  Eq.  (4)  with  a(t)  =  AsinQr.  Figure  6  shows  some  typical  periodic  responses  from  the  system  in  Fig.  5  (numerical 
simulation  in  parts  (a)  and  (b)  and  experimental  in  parts  (c)  and  (d)).  Each  time  series  was  produced  with  a  forcing  amplitude 
of  A  =  6.25  cm  (compare  with  the  Asnap  ss  8.21  cm),  however  parts  (a)  and  (c)  corresponded  to  a  forcing  frequency  of 
Q  =  4.40  rad/s  while  parts  (b)  and  (d)  corresponded  to  Q  3.36  rad/s.  The  transmissibility  of  the  forcing  mechanism 
effectively  translates  to  a  moving  restoring  force,  but  superimposing  the  phase  trajectories  in  this  way  gives  a  clear  indication 
of  whether  the  system  is  oscillating  about,  or  snapping  between,  the  two  sides  of  the  system.  It  is  interesting  that  there  is 
better  agreement  between  the  numerical  and  experimental  results  for  the  large  snap-through  motion  in  Fig.  6(b)  and  (d)  than 
for  the  small  non  snap-through  motion  in  (a)  and  (c).  It  is  likely  that  this  is  due  to  some  nonlinear  damping  effects  that  were 
not  included  in  the  model,  particularly  Coulomb  damping  (already  discussed  earlier),  which  has  a  more  pronounced  effect 
on  slower  motions. 

Figure  7  shows  several  time  series  from  simulation  (a-d)  and  experiment  (e-h),  this  time  for  a  more  shallow  angle 
0o  =  15.0°  (0i  =  29.9°)  and  A  =  1.10  cm  (compare  with  Asltap  k,  1.56  cm).  The  relatively  simple  ’period-one  (PI)’,  periodic 
responses  shown  in  parts  (a,b,e,f)  correspond  to  a  forcing  frequency  of  Q  =  4.9  rad/s,  i.e.,  they  co-exist  at  this  set  of  parameter 
values,  and  which  attractor  occurs  depends  on  the  initial  conditions.  This  aspect  of  the  behavior  will  be  revisited  a  little  later. 

The  non-simple  but  still  periodic  responses  shown  in  parts  (c  and  g)  are  not  quite  the  same.  The  numerical  result  shown 


in  part  (c)  occurred  for  £2  =  7.6  rad/s  and  can  be  characterized  as  a  ’period-four  (P4)’  oscillation,  i.e.,  it  repeats  itself  every 
four  forcing  cycles.  The  experimental  result  shown  in  part  (g)  is  for  a  slightly  different  forcing  frequency  (£2  =  7.9  rad/s) 
and  is  in  fact  a  P5  oscillation.  The  chaotic  behavior  for  the  numerical  result  (part  d)  corresponds  to  £2  -  7.8  rad/s  and  the 
experimental  result  (part  h)  corresponds  to  £2  =  8.1  rad/s.  These  values  are  very  close  to  those  used  to  generate  the  time 
series  in  parts  (c  and  g),  thus  suggesting  a  parameter  sensitivity  in  this  range. 

The  frequency  content  of  these  time  series  (see  Fig.  8  for  the  more  interesting  time  series)  can  be  quite  useful  in 
assessing  certain  qualitative  features  of  the  dynamics.  They  can  also  assist  in  determining  whether  the  system  is  exhibiting 
snap-through  oscillations,  since  it  is  typically  somewhat  difficult  to  consistently  distinguish  between  trajectories  that  snap- 
through  from  those  that  do  not.  The  broadband  nature  of  the  discrete  Fourier  transform  (DFT)  of  a  time  series  is  of  course  a 
well-established  signature  of  chaos. 

Figure  9  contains  a  ’bifurcation’  diagram  of  snap-through  behavior  as  a  function  of  forcing  frequency  with  (0o  =  15.0°, 
A  =  1.10  cm) .  Parts  (a)  and  (b)  show  the  number  of  snap-through  events  per  forcing  cycle  from  simulation  and  experimental 
results,  respectively.  A  snap-through  event  was  considered  to  be  any  crossing  of  the  state  0  =  0.  This  definition  of  snap- 
through  works  relatively  well  for  the  system  because  the  unstable  equilibria  does  not  change  significantly  with  the  forcing. 
Each  data  point  in  part  (a)  indicates  the  result  of  one  simulation  with  the  snap-through  events  per  forcing  cycle  being 
averaged  over  1000  forcing  cycles  after  the  transients  were  given  sufficient  time  to  decay.  The  plot  was  created  by  running 
10  simulations  with  random  initial  conditions  for  100  forcing  frequencies  to  capture  co-existing  attractors.  The  experimental 
data  points  were  instead  obtained  by  a  performing  a  moving  average  of  the  number  snap-through  events  per  forcing  cycle  on 
a  sweep-up  followed  by  a  sweep-down  through  forcing  frequency.  The  sweep-up  and  -down  primarily  differ  in  the  region 
below  £2  =  6.2  rad/s  where  the  sweep-up  did  not  snap-through  and  the  sweep-down  did.  The  responses  that  snap-through 
twice  per  forcing  cycle  indicate  a  snap-through  in-phase  with  the  forcing,  i.e.,  PI  snap-through  (see  Fig.  7(b  and  f)  for 
example  time  series).  The  responses  with  intermediate  occurrences  of  snap-through  indicate  responses  with  occasional,  or 
sporadic,  snap-through,  which  may  occur  in  higher-period  periodic  behavior  (Fig.  7(c  and  g))  and  chaotic  behavior  (Fig.  7(d 
and  h)).  Figure  9(c)  shows  the  likelihood  (via  simulation  only)  of  the  system  response  converging  to  one  of  the  three 
competing  types  of  behavior  seen  it  parts  (a)  and  (b),  where  green  corresponds  to  no  snap-through,  red  corresponds  to 
PI  snap-through,  and  blue  is  higher-periodic  or  chaotic  (less  frequent)  snap-through.  The  likelihood  was  measured  by 
running  100  simulations  at  random  initial  conditions  at  each  frequency.  This  method  may  be  more  difficult  to  apply  to 
higher  dimensional  systems  as  the  separatrices  are  typically  more  complex  and  may  be  more  dependent  on  the  system 
forcing.  However,  in  this  relatively  simple  system  it  gives  a  straightforward  interpretation  of  a  potentially  onerous  event 
(snap-through)  and  thus  useful  data  to  inform  practical  issues  including  life-time  fatigue  predictions. 

Another  approach  might  be  to  observe  some  kind  of  energy  measure  of  a  given  response.  For  a  system  to  snap-through 
it  is  postulated  that  the  external  forcing  must  add  a  sufficient  amount  of  energy  to  the  system.  Figure  10  shows  the  average 
kinetic  energy  over  a  span  of  forcing  cycles  against  the  forcing  frequency  for  the  same  parameters  as  Fig.  9.  The  experimental 
data  points  were  obtained  by  first  sweeping  slowly  up,  then  down,  through  forcing  frequency.  The  agreement  between  the 
experimental  and  numerical  results  is  excellent,  with  only  a  slight  horizontal  shift  in  the  data  points. 


One  anticipates  that  the  average  kinetic  energy  of  a  system  would  be  much  larger  for  trajectories  that  traverse  across 
the  unstable  equilibria  than  for  those  that  oscillate  about  a  single  equilibria;  however  this  detection  is  made  more  difficult 
by  responses  that  only  occasionally  snap-through.  This  typically  occurs  with  very  high  periodic,  or  more  often,  chaotic 
responses,  which  prevail  for  much  of  the  right  half  of  the  energy  plots.  This  makes  the  determination  of  chaos  another 
important  aspect  of  the  analysis,  especially  since  it  is  known  that  single  well  chaos  is  relatively  rare  [24], 

The  standard  method  of  identifying  chaos  is  to  determine  the  Lyapunov  spectrum  of  the  system.  Lyapunov  exponents 
(LE)  are  however  relatively  difficult  to  determine  from  experimental  data  [25] .  Figure  1 1(a)  shows  the  LE  obtained  from  the 
simulation  of  equation  (4)  (using  the  method  in  [26])  against  forcing  frequency.  The  data  points  were  produced  by  running 
10  simulations  for  different  random  initial  conditions  at  100  frequencies  to  make  it  possible  to  capture  co-existing  attractors. 
When  comparing  with  Fig.  10  it  appears  that  most  of  the  chaotic  responses  of  the  system  occur  for  snap-through  trajectories. 

An  alternative  approach  is  introduced  in  [27],  whereby  the  number  of  distinct  frequency  spikes  (peaks  on  the  DFT, 
above  a  certain  threshold)  of  a  dynamical  system  are  plotted  against  the  parameter  used  for  the  simulation/experiment.  There 
is  typically  a  clear  dichotomy  in  the  number  of  peaks  for  a  chaotic  response  with  a  broadband  spectrum  and  a  non-chaotic 
response  with  a  finite  number  of  peaks,  making  identification  of  chaos  possible.  Figure  11(b)  shows  the  identification  of 
chaos  using  this  peak  count  method.  In  these  results  10  simulations  with  randomly  generated  initial  conditions  were  run 
at  100  forcing  frequencies.  It  is  confirmed  that  the  chaotic  responses  are  almost  entirely  within  the  region  of  high  energy 
response. 

An  interesting  ’second-order’  feature  in  these  plots  is  an  assessment  of  how  chaotic  a  signal  might  be.  Two  sample  time 
series  (indicated  by  the  vertical  dashed  red  lines)  are  also  shown  in  Fig.  11(c)  and  (d).  Part  (d)  corresponding  to  £2  =  7.6 
rad/s  shows  a  slightly  stronger  periodicity  than  the  time  series  shown  in  part  (c)  (£2  =  7.0  rad/s):  the  second  half  of  the 
time  series  appears  nearly  periodic,  this  short  term  periodicity  is  ubiquitous  throughout  the  time  series  as  the  trajectory 
appears  to  spend  more  time  oscillating  about  (rather  than  between)  the  underlying  static  equilibria.  When  subject  to  the 
LE  algorithm  this  periodic  signature  results  in  a  slightly  less  positive  LE,  and  there  are  a  few  less  spectral  peaks  occurring 
above  the  threshold.  It  is  not  uncommon  for  a  chaotic  trajectory  to  spend  periods  of  time  in  a  near-periodic  state  (a  typical 
chaotic  attractor  has  embedded  within  it  many  unstable  periodic  orbits).  It  was  also  just  mentioned  that  in  certain  (relatively 
small)  regimes  of  forcing  frequency  co-existing  oscillations  occur.  In  order  to  illuminate  this  feature.  Fig.  11(e)  shows  the 
outcome  of  100  simulations  using  randomly  generated  initial  conditions,  in  which  the  percentage  of  numerical  simulations 
that  lead  to  chaos  is  charted  for  both  the  LE  (dashed  black  line)  and  the  peak-count  (solid  green  line)  approaches.  Clearly, 
both  criteria  are  effective  in  distinguishing  chaotic  from  non-chaotic  motion.  Outside  of  this  range  of  forcing  frequencies  the 
behavior  was  typically  non-chaotic,  although  co-existing  attractors  could  still  exist,  and  the  periodic  motion  could  be  either 
small-amplitude  (contained  within  the  vicinity  of  an  equilibrium),  or  relatively  large-amplitude  (cross- well)  motion. 

The  same  approach  was  also  applied  to  the  experimental  system.  Figure  12  shows  bifurcation  diagrams,  based  on  18 
tests,  in  terms  of  chaotic  vs.  non-chaotic  behavior.  The  computation  of  the  largest  Lyapunov  exponent  (A.)  is  based  on  a 
standard  algorithm  [28]  in  which  a  time  series  is  assessed  in  terms  of  the  exponential  rate  at  which  nearest  neighbors  tend 
to  evolve  in  time.  This  can  be  a  sensitive  undertaking,  for  example,  part  (a),  for  £2  =  6.518,  and  (b),  for  Q  =  8.140  of 


Fig.  12  show  how  the  approach  requires  a  linear  fit.  This  is  by  no  means  straightforward.  The  red  dashed  lines  in  (a)  and  (b) 
show  an  approximately  exponential  divergence  (the  y-axis  is  a  log-scale)  for  short  time  evolutions.  However,  the  response 
in  part  (a)  is  characterized  by  a  very  small  slope  and  corresponds  to  a  periodic  oscillation,  whereas  the  case  in  part  (b)  is 
more  convincingly  positive  and  corresponds  to  a  chaotic  oscillation.  Part  (c)  of  this  figure  summarizes  this  behavior,  and 
part  (d)  shows  the  equivalent  results  based  on  the  peak  count  criterion.  The  red  dashed  lines  in  these  parts  show  the  specific 
frequencies  at  which  the  snapshots  were  taken  for  parts  (a)  and  (b).  The  peak  counts  approach  appears  to  be  more  robust  and 
easier  to  apply  than  the  conventional  LE  approach.  The  agreement  with  the  numerical  results  is  quite  good.  The  bifurcations 
are  slightly  shifted,  however  the  progression  and  even  some  of  the  periodic  windows  are  still  captured. 


6  Conclusions 

Snap-through  buckling  is  investigated  on  a  SDOF  system.  Some  aspects  of  the  underlying  force-deflection  and  free 
vibration  behavior  of  the  system  are  described.  Under  the  action  of  periodic  excitation,  the  focus  of  interest  shifts  from 
equilibria  to  oscillatory  behavior.  In  a  broad  sense,  the  behavior  of  the  forced  system  can  be  divided  into  regimes  of  either 
small-amplitude  or  large-amplitude  behavior.  This  latter  response  can  be  associated  with  snap-through  behavior  -  an  often 
highly  undesirable  behavior  in  many  structural  applications.  Alternatively,  the  behavior  can  be  divided  into  regimes  of 
chaotic  and  non-chaotic  behavior.  Chaos  is  commonly  associated  with  snap-through  behavior  but  they  are  not  necessarily 
the  same  thing.  A  promising  method  of  distinguishing  between  small-  and  large-amplitude  responses  based  on  average  total 
energy  is  introduced  and  applied  to  the  system  with  good  agreement  between  experimental  results  and  the  numerical  model. 
Another  method  based  on  distinguishing  between  chaotic  and  non-chaotic  responses  is  also  proposed.  In  general,  this  work 
provide  a  useful  overview  of  how  snap-through  behavior  depends  on  the  forcing  parameters. 
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Figure  Captions: 


1 .  A  typical  curved  panel,  and  two  schematic  scenarios  in  which  such  a  system  might  exhibit  a  snap-through  event  in  its 
force-displacement  relationship,  (a)  pre-snap,  (b)  post-snap,  (c)  limit  point  buckling,  (d)  pitchfork  bifurcation. 

2.  A  single  degree  of  freedom  (SDOF)  link  model. 

3.  (a)  Photograph  of  experimental  setup,  (b)  the  low-friction  pin  joint,  (c)  the  Scotch-yoke  forcing  mechanism. 

4.  Identification  of  damping  parameter  |3  (Kg/s),  (a)  A  typical  nonlinear  free  decay;  the  points  are  experimental  data, 
and  the  continuous  line  represents  the  numerical  integration  of  equation  (4)  with  (1  =  1.2.  (b)  Normalized  average  error  vs. 
p  for  the  large  amplitude  time  series  (in  part  (a)),  (c)  Normalized  average  error  vs.  (1  for  a  small  amplitude  time  series. 

5.  Free  response  characteristics,  (a)  force  vs.  natural  frequency  (squared),  (b)  force  vs.  deflection,  (c)  natural  frequency 
(squared)  vs.  deflection.  The  points  are  experimental  data,  the  continuous  lines  are  the  theoretical  results. 

6.  Experimental  and  simulated  time  series  superimposed  on  the  restoring  force.  Numerical,  (a)  Q  =  4.40  rad/s;  (b) 
Q  =  3.36  rad/s;  Experimental,  (c)  Q  =  4.40  rad/s;  (d)  Q  =  3.36  rad/s. 

7.  Numerical  (a-d)  and  experimental  (e-h)  time  series,  (a)  Q  =  4.9  rad/s;  (b)  Q  =  4.9  rad/s;  (c)  Q  =  7.6  rad/s;  (d) 
Q  =  7.8  rad/s;  (e)  Q  =  4.9  rad/s;  (f)  Q  =  4.9  rad/s;  (g)  Q  =  7.9  rad/s;  (h)  Q  =  8.1  rad/s. 

8.  Experimental  and  simulated  DFT’s  for  parts  (c,d,g,h)  in  Fig.7  respectively. 

9.  Occurrence  of  snap-through,  (a)  simulation,  (b)  experiment,  (c)  relative  dominance  of  co-existing  attractors  (simu¬ 
lation  only).  Green  -  non-snap,  red  -  PI  snap-through,  blue  -  higher  periodic  or  chaotic  (less  frequent)  snap-through.  The 
vertical  dashed  red  lines  in  (a)  and  (b)  indicate  the  specific  frequencies  relating  to  Fig.7. 

10.  Average  kinetic  energy  as  a  function  of  forcing  frequency,  (a)  simulation,  (b)  experiment. 

11.  Distinction  between  chaotic  and  non-chaotic  behavior  based  on  (a)  the  largest  Lyapunov  exponent,  (b)  the  peak- 
count  criterion;  (c)  and  (d)  typical  chaotic  time  series,  (e)  relative  dominance  of  chaotic  behavior. 

12.  Experimental  LE  and  peak  count,  (a)  and  (b)  typical  linear  fits  for  the  local  rate  of  divergence,  (c)  Largest  LE  as  a 
function  of  the  forcing  frequency,  (d)  corresponding  peak  count  result. 
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A  Heuristic  Method  for  Identifying  Chaos  from  Frequency  Content 

R.  Wiebe1’11)  and  L.N.  Virgin1’15) 

School  of  Engineering,  Duke  University,  Durham,  NC. 
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The  sign  of  the  largest  Lyapunov  exponent  is  the  fundamental  indicator  of  chaos  in  a  dynamical  system. 
However,  although  the  extraction  of  Lyapunov  exponents  can  be  accomplished  with  (necessarily  noisy)  exper¬ 
imental  data,  this  is  still  a  relatively  data-intensive  and  sensitive  endeavor.  This  paper  presents  an  alternative 
pragmatic  approach  to  identifying  chaos  using  response  frequency  characteristics  and  extending  the  concept 
of  the  spectrogram.  The  method  is  shown  to  work  well  on  both  experimental  and  numerical  data. 
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It  is  often  easy  to  distinguish  between  chaotic 
and  nonchaotic  trajectories  by  visual  inspection 
of  their  time  series  and  frequency  spectra.  The 
frequency  transform  of  chaotic  time  series  gener¬ 
ally  show  certain  characteristic  signatures  such  as 
fractal  geometry  and  broadband  frequency  con¬ 
tent.  This  is  contrasted  with  the  distinct  sig¬ 
nature  of  periodic,  or  nonchaotic  behavior.  The 
heuristic  method  developed  herein  uses  a  quan¬ 
tifiable  measure  based  on  the  broadband  non¬ 
smooth  nature  of  the  frequency  spectra  of  chaotic 
trajectories.  The  method  is  then  verified  numeri¬ 
cally  and  experimentally  by  comparison  with  the 
largest  Lyapunov  exponent. 


I.  BACKGROUND 

Chaos  is  defined  as  the  output  of  a  deterministic  dy¬ 
namical  system  that  exhibits  sensitivity  to  initial  con¬ 
ditions.  If  the  governing  equations  are  known  then  it  is 
usually  straightforward  to  interrogate  (especially  chaotic) 
behavior.  In  the  experimental  (noisy)  context,  however, 
it  is  more  challenging  to  characterize  behavior  under  cer¬ 
tain  parameter  regimes.  An  excellent  summary  of  time 
series  analysis  and  chaos  is  provided  in1-3.  The  Lya¬ 
punov  exponent  (LE)  spectrum  of  a  dynamical  system 
describes  the  divergence/convergence  rate  of  nearby  tra¬ 
jectories  in  state  space4  and  is  the  most  direct  measure 
of  chaos.  Determining  the  LE  spectrum  of  a  numeri¬ 
cal  system  can,  in  principle,  be  done  directly  as  the  full 
system  state  is  known.  An  algorithm  for  calculating  the 
LE’s  of  a  numerical  system  is  given  in5.  For  experimental 
data,  however,  the  process  is  more  difficult  since  not  all 
states  are  typically  measured  and  because  the  lineariza¬ 
tion  (or  normalization)  required  in  calculating  the  LE  is 
not  possible.  Therefore  in  order  to  extract  the  LE’s  of  an 
experimental  system  one  must  normally  reconstruct  the 
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phase  portrait  via  time  lag  embedding  and  then  use  the 
method  of  nearest  neighbors  to  analyze  the  divergence 
rates6-8.  This  method  typically  requires  long  time  series 
to  ensure  accuracy  and  requires  some  tuning  of  the  time 
span  over  which  the  divergence  remains  approximately 
exponential,  i.e.  the  time  span  over  which  the  LE  may 
be  calculated.  The  difficulty  of  calculating  LE’s  from 
noisy  time  series  is  discussed  in9. 

In  many  dynamical  systems  the  presence  of  noise  can 
lead  to  a  phenomenon  known  as  unstable  dimension  vari¬ 
ability  (see10  for  an  overview).  This  occurs  when  the 
number  of  stable,  unstable,  and  neutrally  stable  eigendi- 
rections  change  intermittently  when  a  system  is  per¬ 
turbed  between  co-existing  invariant  sets  by  noise.  One 
anticipates  that  this  phenomenon  would  further  compli¬ 
cate  the  calculation  of  LE’s. 

Another  interesting  approach  for  identifying  chaos  is 
discussed  in11  whereby  noise  is  artificially  added  to  a  sig¬ 
nal  in  order  to  find  the  threshold  noise  limit  whereby 
the  nonlinearity  becomes  undetectable.  This  ’titration’ 
method  is  shown  to  work  quite  well  and  has  the  added 
benefit  of  also  revealing  the  level  of  noise  in  the  system, 
however  it  also  requires  one  to  choose  some  nonlinearity 
indicator. 

Frequency  analysis  is  commonly  used  to  characterize 
noisy  signals  and  is  a  standard  procedure  in  signal  pro¬ 
cessing.  The  power  spectrum  of  noisy  signals  has  long 
been  used  to  characterize  the  noise  type  with  the  power 
spectrum  power  law12.  The  frequency  content  of  deter¬ 
ministic  or  systems  with  minimal  noise  is  also  commonly 
used  to  qualitatively  identify  whether  a  time  series  is 
chaotic  or  periodic13-15.  However,  the  inspection  must 
be  done  qualitatively  and  therefore  the  method  is  not 
tractable  for  a  full  investigation  of  the  parameter  space 
of  a  dynamical  system.  Attempts  have  been  made  to 
automate  the  identification  process  by  assigning  a  single 
metric  to  the  power  spectrum  or  DFT  of  a  time  history 
such  as  its  fractal  dimension16.  The  heuristic  method  dis¬ 
cussed  and  developed  herein,  which  consists  of  a  counting 
the  peaks  on  the  DFT,  was  first  mentioned  in1'.  How¬ 
ever,  the  method  was  only  applied  to  a  numerical  system, 
and  the  effects  of  noise  were  not  thoroughly  discussed. 
The  effect  of  noise  on  the  method  was  studied  numeri¬ 
cally  in18.  In  the  following,  the  same  method  is  applied 
to  an  experimental  system,  which  serves  to  validate  its 
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practicality,  and  confirm  its  ability  to  work  well  even  in 
the  presence  of  real  noise. 


II.  THE  HEURISTIC  METHOD 

It  is  often  straightforward  to  distinguish  between 
chaotic  and  nonchaotic  trajectories  by  inspecting  their 
frequency  spectra.  A  periodic  signal  will  have  a  spec¬ 
trum  dominated  by  a  single  peak  (harmonic)  or  a  finite 
set  of  peaks  for  a  more  complicated  but  still  periodic, 
or  quasi-periodic,  signal.  Noise  tends  to  cause  a  slight 
spread  in  the  width  of  these  peaks  and  some  underly¬ 
ing  low  level  energy  across  the  frequency  spectra,  but 
the  peaks  are  generally  still  easy  to  distinguish.  For  a 
chaotic  response,  on  the  other  hand,  the  spectral  content 
is  spread  over  a  relatively  broad  range  of  frequencies.  A 
spectrogram  or  waterfall  plot19  summarizes  these  effects 
as  a  system  parameter  is  changed,  i.e. ,  there  is  a  distinct 
change  in  the  frequency  signature  at  a  bifurcation.  In 
some  circumstances  we  wish  to  display  different  types  of 
behavior  as  two  parameters  are  changed.  Spectrograms 
are  however  limited  to  showing  bifurcations  in  a  single 
parameter,  as  two  dimensions  are  needed  to  display  the 
frequency  spectra  themselves. 

The  Heuristic  method  used  herein  extends  the  idea  of 
the  spectrogram  to  allow  plotting  of  multi-parameter,  es¬ 
pecially  two-parameter,  systems.  The  heuristic  approach 
relies  on  the  noisy  nature  of  the  Fourier  transform  of 
chaotic  trajectories  and  is  summarized  in  Fig.  1.  First 
a  DFT/FFT  is  applied  to  a  time  series  of  the  dynami¬ 
cal  system  of  interest.  For  a  steady-state  response  it  is 
appropriate  to  use  just  a  single  sample  (as  used  in  the  ex¬ 
amples  in  this  paper) ,  providing  sufficient  data  is  used  to 
capture  the  essential  dynamics.  Just  as  in  the  calculation 
of  LE’s  one  needs  to  allow  the  system  to  settle  onto  the 
attractor  prior  to  beginning  the  calculation  of  the  DFT. 
Once  the  DFT  is  obtained  a  threshold  value  is  chosen 
based  on  the  maximum  peak  height  of  the  frequency  spec¬ 
tra,  represented  by  the  dashed  line  in  Fig.  1.  The  number 
of  peaks  (local  maxima)  above  this  threshold  are  counted. 
If  the  number  of  peaks  counted  is  above  a  second  thresh¬ 
old  the  system  is  considered  chaotic,  otherwise  it  is  la¬ 
beled  periodic.  The  challenge  arises  in  choosing  the  two 
threshold  values.  The  DFT  of  a  periodic  signal  contains 
only  a  handful  of  peaks  at  resonant  frequencies  as  in  Fig. 
1(a).  However  in  the  presence  of  noise  (or  rounding  error 
in  the  case  of  simulation)  the  curve,  when  magnified,  will 
contain  many  local  maxima  especially  near  zero  (as  in 
Fig.  1(b)).  The  first  threshold  can  therefore  usually  be 
chosen  to  be  a  small  percentage  of  the  maximum  peak 
height,  and  therefore  only  the  excited  frequencies  will  be 
counted  for  a  periodic  signal.  Chaotic  signals,  however, 
produce  many  peaks  on  a  similar  order  of  magnitude  to 
the  maximum  peak  as  seen  in  Fig.  1(c).  The  dichotomy 
of  peak  counts  seen  in  Fig.  1  are  typical  of  many  systems 
and  usually  make  it  relatively  easy  to  tune  the  threshold 
values. 
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FIG.  1.  Illustration  of  heuristic  method  with  typical  signals, 
(a)  periodic,  (b)  periodic  with  noise,  (c)  chaotic. 


Typically  any  state  of  a  system  may  be  used  to  gen¬ 
erate  the  DFT  as  periodic  motion  in  one  state  of  a  cou¬ 
pled  system  typically  implies  periodic  motion  in  the  other 
states.  However,  certain  states  may  produce  frequency 
spectra  in  which  chaos  is  easier  to  ascertain,  for  instance 
in  a  pendulum  angular  velocity  is  bounded  whereas  the 
angle  is  not. 

In  essence  this  method  emulates  visual  identification 
of  chaos  from  frequency  spectra,  and  therefore  allows  for 
high  resolution  sweeps  through  system  parameters  in  a 
relatively  hands-off  fashion.  This  allows  one  to  build  an 
indicator  function  over  the  parameter  space  with  a  bi¬ 
nary  output  describing  the  response  as  either  chaotic  or 
periodic.  As  will  be  seen,  no  distinction  is  made  between 
different  types  of  periodic  behavior.  In  similarity  with 
any  bifurcation  diagram,  capturing  coexisting  attractors 
would  require  additional  testing  based  on  varying  initial 
conditions.  In  practice  however,  the  method  is  most  use¬ 
ful  for  plotting  boundaries  in  two-parameter  space  sepa¬ 
rating  chaotic  from  nonchaotic  motion. 
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One  particular  advantage  of  this  method  over  LE’s  is 
the  ease  of  implementation;  one  needs  only  apply  a  fre¬ 
quency  transform,  a  standard  signal  processing  feature. 
It  may  also  be  used  together  with  LE’s  to  provide  a  level 
of  mutual  verification.  Despite  its  simplicity  the  method 
proves  to  be  surprisingly  robust  (in  the  presence  of  real 
noise)  and  matches  excellently  with  the  LE’s  of  the  nu¬ 
merical  model,  while  being  much  easier  to  obtain  than 
experimental  LE’s. 

The  method  may  run  in  complications  for  systems  with 
a  low  signal  to  noise  ratios  due  to  the  noise  raising  the 
number  of  peaks  on  the  DFT.  One  possible  method  to 
counteract  this  is  discussed  in1,  where  several  DFT’s  are 
averaged  to  reduce  noise  in  the  DFT. 


III.  TEST  CASE:  A  NONLINEAR  MECHANICAL 
OSCILLATOR 

A.  Experimental  Setup 

The  heuristic  method  was  tested  on  the  experimental 
nonlinear  mechanical  oscillator  shown  in  Fig.  2.  The 
oscillator  is  a  product  of  ’’Pasco  Systems”  and  consists 
of  a  pendulum  with  a  torque  spring.  The  pendulum  is 
comprised  of  an  aluminum  disk,  with  an  off-center  brass 
weight.  The  torque  spring  is  made  up  of  a  pair  of  lin¬ 
ear  springs  connected  by  a  string  that  is  looped  over  a 
pulley  attached  behind  the  pendulum.  The  springs  also 
supply  the  forcing  to  the  system.  The  forcing  is  trans¬ 
mitted  by  the  active  spring  attached  to  the  pendulum 
by  the  displacement  of  the  forcing  arm  on  the  motor,  as 
shown  in  the  bottom  inset  of  Fig.  2.  The  forcing  motor 
is  connected  to  the  active  spring  by  a  string  which  passes 
through  the  guide-post  (also  visible  in  the  bottom  inset 
of  Fig.  2).  The  guide-post  is  necessary  to  transform  the 
circular  motion  of  the  forcing  arm  into  a  linear,  nearly 
sinusoidal,  motion  on  the  active  spring.  The  secondary, 
anchor,  spring  ensures  that  both  springs  are  always  in 
tension  to  avoid  needing  a  tension-compression  capable 
spring.  The  response  of  the  system  is  measured  by  a 
rotary  motion  sensor  attached  to  the  axle  of  the  pendu¬ 
lum.  The  forcing  is  also  measured  by  a  rotary  motion 
sensor  attached  to  the  forcing  string  directly  above  the 
guide-post. 

The  potential  energy  function  underlying  the  system  is 
the  summation  of  a  sinusoidal  gravity  potential  due  to  the 
eccentric  mass  (the  brass  attachment),  and  a  parabolic 
potential  due  to  the  linear  torque  spring.  By  altering  the 
eccentric  mass  and  the  torque  spring  stiffness  one  is  able 
to  choose  the  number  of  static  equilibria  of  the  system; 
one  may  choose  1, 3,  5, ...  equilibira  in  an  appropriate  op¬ 
erational  range.  Another  parameter  in  the  experimental 
system  is  the  zero  point  of  the  forcing,  i.e.  the  angle 
of  the  forcing  arm  and  the  corresponding  position  of  the 
static  equilibria.  For  this  paper  the  spring  stiffness  and 
mass  were  chosen  such  that  there  were  three  equilibria 
with  the  forcing  zeroed,  with  the  unstable  equilibria  be- 


FIG.  2.  Experimental  apparatus. 


ing  the  inverted  pendulum  state  and  therefore  two  sym¬ 
metric  stable  equilibria  (one  of  which  is  shown  in  the  top 
inset  figure  2).  The  force  zero  point  was  chosen  to  be 
when  the  forcing  arm  was  horizontal.  This  was  selected 
because  setting  a  horizontal  forcing  arm  coinciding  with 
a  vertical  unstable  equilibrium  (i.e.  inverted  pendulum 
state)  was  the  easiest  to  consistently  calibrate  between 
experimental  runs. 


B.  Numerical  Modelling 

Figure  3  shows  a  schematic  of  the  experimental  appa¬ 
ratus,  where  9  is  the  angle  of  a  ray  projected  to  the  ec¬ 
centric  mass  from  the  pivot  measured  counter-clockwise 
from  vertical,  yo(t)  is  the  time  varying  linear  displace¬ 
ment  of  the  driving  string  above  the  guide-post,  m  is 
the  mass  of  the  brass  attachment,  M  is  the  mass  of  the 
aluminum  disc,  rm  is  the  radius  of  the  aluminum  disc 
(note  that  the  eccentric  mass  is  centered  at  the  edge  of 
the  disc),  rp  is  the  radius  of  the  pulley  attaching  the 
springs  to  the  pendulum,  k\  and  &2  are  the  stiffnesses  of 
the  active  and  anchor  springs  respectively,  and  d  is  the 
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distance  between  the  guide-post  and  the  axle  of  the  driv¬ 
ing  motor.  The  two  adjustable  forcing  parameters  for 
the  system  are  the  length  of  the  forcing  arm  L  and  the 
forcing  frequency  /.  The  rotary  motion  sensor  attached 
to  the  forcing  string  (located  between  the  guide-post  and 
spring  1)  is  omitted  for  clarity.  The  equation  of  motion 
describing  the  system  is  given  by 

d  +  2 (uq9  +  ulO  -  mg^m  sin#  =  -y^y0(t), 

=  (fci  +k2)rp2,  (1) 

I  =  mrm 2  +  ^ Mrm 2, 

where  a  small  linear  viscous  damping  term  with  coeffi¬ 
cient  (  is  inserted  to  capture  the  viscous  energy  dissipa¬ 
tion  of  the  system. 


FIG.  3.  Schematic  of  system  (not  to  scale). 

The  choice  of  zero-ing  the  forcing  with  the  forcing  arm 
horizontal  leads  to  a  slightly  biased  forcing  due  to  the 
geometry  relating  the  forcing  arm  and  the  guide-post  as 
shown  in  Fig.  3.  The  motion  of  the  string  above  the 
guide-post  driving  the  active  spring  is  given  by 


yo(t)  =  \/d2  —  2  Ld  cos  (2-Kft)  +  L2  —  ^  d2  +  L2.  (2) 

Note  that  it  was  decided  to  not  normalize  the  length 
scale  as  the  units  are  of  a  physically  appreciable  scale, 
and  there  is  no  clear  meaningful  characteristic  length  for 
the  system.  Figure  4  shows  a  time  series  and  a  time 
lag  plot  (with  a  time  lag  of  one  forcing  period)  of  Eqn. 
(2)  versus  experimental  measurement  and  a  sinusoid  for 
comparison.  It  is  clear  that  Eqn.  (2)  captures  both  the 
bias  and  the  shape  better  than  the  sinusoid. 


yo(t)(cm)  (a) 


y0(t+T)(cm)  (b) 


FIG.  4.  System  identification:  forcing.  Dots  =  experimental 
data,  black  line  =  Eqn.  (2),  red  line  =  sinusoid,  (a)  Time 
series,  (b)  time  lag  portrait. 

Damping  is  provided  by  eddy  current  drag  from  the 
interaction  between  the  aluminum  disc  and  a  permanent 
magnet  placed  slightly  behind  it  (visible  in  the  top  inset 
of  Fig.  2).  This  form  of  damping  is  described  in20  and  is 
shown  to  be  linear  viscous  damping  as  long  as  velocities 
are  not  too  large.  Determining  at  which  velocity  the  lin¬ 
ear  viscous  damping  assumption  fails  requires  extensive 
electromagnetic  analysis  and  was  not  investigated.  The 
excellent  agreement  between  experimental  and  numerical 
results  shown  later,  however,  suggest  that  the  velocities 
obtained  experimentally  are  well  below  this  threshold. 

For  a  linear  system  there  are  many  well-tested  meth¬ 
ods  of  assessing  energy  dissipation  such  as  logarithmic 
decrement  and  the  half-power  methods21.  In  order  to  ap¬ 
ply  these  methods  the  eccentric  brass  mass  was  removed 
leaving  a  simple  linear  oscillator  (Eqn.  (1)  with  m  =  0). 
The  damping  coefficient  obtained  by  removing  the  mass 
must  however  be  corrected  by 

C  =  \j ^pCm=o  »  0.90Cm=o.  (3) 

The  correction  is  necessary  because  the  damping  is  lin¬ 
early  proportional  to  velocity  prior  to  normalization,  i.e. 
the  term  2^ojqI  is  constant.  Therefore  £  must  change 
inversely  with  \fl  since  u> o  oc  1  /y/l. 

Figure  5  shows  the  results  of  three  different  methods  of 
assessing  damping.  The  logarithmic  decrement  method 
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uses  successive  peak  heights  in  the  free  decay  time  series’ 
as  shown  in  5  (a).  The  damping  coefficient  can  be  ap¬ 
proximated  by  the  relation  £  =  In  (An / An+m) /27tto.  The 
blue  and  black  lines  in  part  (a)  both  yielded  Cm=o  ~  6.0% 
or  £  ss  5.4%  for  the  corrected  damping  coefficient.  Part 
(b)  shows  the  frequency  response  curve  for  a  forcing  am¬ 
plitude  of  L  =  5.5 cm.  This  curve  is  used  in  the  half 
power  method  to  obtain  the  damping  coefficient  from 
the  relation  £  =  (/2  —  /i)/2/0  which  yields  £m=o  ~  6.5% 
or  £  «  5.9%  for  this  system.  A  subtle  issue  with  this 
approach  is  that  the  half  power  relation  is  derived  for 
a  sinusoidal  forcing  function,  however  as  already  shown 
the  forcing  for  this  system  is  not  quite  sinusoidal.  Fi¬ 
nally  part  (c)  shows  the  result  of  an  error  minimization 
approach,  whereby  the  time  averaged  error  between  a 
simulated  free  decay  (Eqn.  (1)  with  m  =  0,  yo  =  0) 
and  the  experimental  time  series  (from  part  (a))  is  plot¬ 
ted  against  the  £m=o  value  used  for  the  simulation.  The 
minimum  error  occurs  for  £m= o  ~  6.2%  or  £  «  5.6%.  The 
first  inclination  when  presented  with  multiple  measure¬ 
ments  is  to  take  their  average.  However,  since  the  three 
methods  are  so  different,  it  is  not  necessarily  true  that 
the  average  will  be  any  more  physically  relevant  than  any 
individual  result.  It  is  also  unlikely  that  the  results  are 
accurate  to  the  tenth  of  a  percent.  Therefore  it  was  de¬ 
cided  to  simply  take  £  =  6%  as  two  of  the  methods  round 
to  this  value. 


for  every  3  forcing  cycles)  oscillation  with  the  addition 
of  the  phase  portraits  to  highlight  the  more  complicated 
motion.  Finally  Fig.  9  shows  the  numerically  and  exper¬ 
imentally  obtained  time  series’,  phase  portraits,  DFT’s, 
and  Poincare  sections  for  a  chaotic  oscillation.  The  phase 
portraits  were  all  obtained  using  time  lag  embedding  with 
a  time  lag  of  a  quarter  of  the  forcing  period.  The  chaotic 
Poincare  section  data  was  however  taken  once  per  forc¬ 
ing  cycle,  and  numerical  differentiation  was  used  to  ob¬ 
tain  the  angular  velocity  experimentally.  (The  dashed 
line  running  across  the  DFT’s  for  all  of  these  samples  is 
at  5%  of  the  maximum  peak  height.  This  is  the  cut-off 
used  for  peak  counting  in  the  Heuristic  method  applied 
later.) 
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FIG.  6.  Small  amplitude  period  1  oscillation,  L  =  5.5 cm  and 
/  =  1.10 Hz.  (a)  Numerical  time  series,  (b)  numerical  DFT, 
(c)  experimental  time  series,  (d)  experimental  DFT. 
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FIG.  5.  System  identification:  damping,  (a)  Logarithmic 
decrement,  (b)  half  power  method  with  /i  =  0.65 Hz,  fo  = 
0.69 Hz,  and  fo  =  0.74 Hz,  (c)  error  minimization. 


The  equation  of  motion  (Eqn.  (1))  was  simulated  using 
a  fourth-order  Runge-Kutta  time  stepping  scheme22,  and 
is  in  excellent  agreement  with  the  experimental  results. 
Figures  6  and  7  show  the  numerical  and  experimentally 
obtained  time  series’  and  DFT’s  of  a  small  amplitude 
and  large  amplitude  period  1  (i.e.  one  response  cycle 
for  every  forcing  cycle)  oscillation  respectively.  Figure  8 
shows  the  same  for  a  period  3  (i.e.  one  response  cycle 


FIG.  7.  Large  amplitude  period  1  oscillation,  L  =  5.5cm  and 
/  =  0.55 Hz.  (a)  Numerical  time  series,  (b)  numerical  DFT, 
(c)  experimental  time  series,  (d)  experimental  DFT. 
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FIG.  8.  Period  3  oscillation,  L  =  5.5cm  and  /  =  0.39Hz.  (a)  Numerical  time  series,  (b)  numerical  phase  portrait,  (c)  numerical 
DFT,  (d)  experimental  time  series,  (e)  experimental  phase  portrait,  (f)  experimental  DFT. 


FIG.  9.  Chaotic  oscillation,  L  =  5.5 cm  and  /  =  0.95 Hz.  (a)  Numerical  time  series,  (b)  numerical  phase  portrait,  (c)  numerical 
DFT,  (d)  numerical  Poincare  section,  (e)  experimental  time  series,  (f)  experimental  phase  portrait,  (g)  experimental  DFT,  (h) 
experimental  Poincare  section. 


C.  Identification  of  Chaos 

There  are  many  applications  where  one  is  interested 
the  parameter  dependence  of  a  system.  There  are  many 
forms  of  bifurcation  diagrams  used  to  present  the  qualita¬ 
tive  changes  (in  our  case  whether  the  response  is  chaotic 
on  nonchaotic)  in  behavior  with  respect  to  a  change  in 
parameters.  Producing  bifurcation  diagrams,  especially 
experimentally,  is  typically  a  laborious  process.  Any 
method  that  would  reduce  the  amount  of  data  required  to 
qualify  the  nature  of  a  response  as  chaotic  or  nonchaotic 
would  therefore  be  highly  advantageous. 

Spectrograms  are  at  the  heart  of  the  heuristic  method 
and  may  be  used  as  a  one-dimensional  bifurcation  di¬ 


agram.  Figure  10  shows  numerical  and  experimental 
spectrograms  over  a  range  of  forcing  frequencies  with  the 
forcing  amplitude  held  constant  at  L  =  5.5 cm.  The  ordi¬ 
nate  of  the  spectrograms  is  the  response  frequency  while 
the  contour  indicates  the  power  of  the  system  (dark  col¬ 
ors  are  low  power,  light  colors  high  power)  at  that  fre¬ 
quency.  The  resolution  of  the  spectrograms  is  100  and  58 
individual  DFT’s  respectively  for  the  numerical  and  ex¬ 
perimental  plots  respectively,  with  each  DFT  containing 
300  data  points  between  0  and  3  Hz.  Spectrograms  are 
useful  in  that  they  clearly  indicate  the  components  of  a 
periodic  response,  and  show  chaotic  regions  as  bands  of 
erratic  contouring.  Their  downside  is  that  since  one  di¬ 
mension  is  used  to  plot  the  response  frequency,  only  one 
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system  parameter  may  be  studied  at  a  time.  The  heuris¬ 
tic  method  works  to  free-up  this  dimension,  rather  than 
plotting  the  entire  DFT  of  a  response,  only  the  number  of 
peaks  is  plotted.  In  doing  this  one  loses  the  information 
about  the  nature  of  a  periodic  response,  however  there 
is  still  a  clear  difference  between  chaotic  and  periodic  re¬ 
sponses.  Figure  11  shows  the  number  of  peaks  out  of  a 
total  possible  300  points  on  the  DFT  versus  the  forcing 
frequency  (again  for  L  =  5.5 cm)  for  the  numerical  and 
experimental  system.  In  order  to  make  the  comparison 
more  direct,  the  numerical  resolution  (again  along  the 
forcing  frequency  axis)  was  reduced  down  to  58  to  match 
the  number  of  experimental  tests  performed.  The  cut-off 
threshold  for  this  plot  was  set  at  1%  of  the  maximum 
peak  height  to  eliminate  false  peaks  (due  to  noise)  being 
counted.  It  is  clear  that  this  method  produces  a  useful 
alternative  bifurcation  diagram  as  the  number  of  peaks  in 
the  chaotic  and  nonchaotic  regions  is  easy  to  distinguish. 
Note  that  dashed  vertical  lines  denote  the  locations  of 
the  four  samples  in  Fig.’s  6  to  9. 

All  of  the  time  series  used,  both  experimentally  and 
numerically,  were  first  allowed  to  run  for  a  period  of  200 
seconds  to  allow  transients  to  decay  and  then  the  follow¬ 
ing  100  seconds  were  recorded  for  analysis.  The  selection 
of  these  two  times  should  be  made  on  a  system  specific 
basis.  In  the  case  of  periodically  forced  systems  this  is 
made  somewhat  easier,  as  one  may  just  select  a  time  pe¬ 
riod  long  enough  to  capture  several  forcing  periods.  How¬ 
ever,  one  also  needs  to  consider  how  many  points  they 
require  on  the  DFT  of  the  response.  This  is  important 
when  applying  the  heuristic  method  because  the  number 
of  peaks  on  a  periodic  response  does  not  typically  increase 
with  resolution,  whereas  for  a  chaotic  response  the  num¬ 
ber  does  increase  with  resolution.  Therefore  there  should 
be  enough  points  to  allow  for  a  clear  difference  in  peak 
counts,  where  this  distinction  will  grow  more  and  more 
clear  with  increasing  resolution. 

The  sensitivity  of  the  peak  count  cut-off  threshold  is 
investigated  in  Fig.  12  for  the  experimental  data.  Part 
(a)  shows  that  a  wide  range  of  threshold  values  from  1% 
to  10%  of  the  maximum  peak  height  yield  easily  distin¬ 
guishable  regions  of  chaos  and  periodic  behavior.  For  the 
case  of  a  0%  cut-off  periodic  responses  begin  to  be  masked 
by  noise,  making  it  difficult  to  distinguish  chaotic  from 
periodic  responses.  This  is  expected  and  is  due  to  the  low 
magnitude  peaks  in  the  DFT  discussed  earlier.  Part  (b) 
shows  the  peak  count  values  versus  the  cut-off  threshold 
for  a  periodic  (/  =  0.55  Hz)  and  chaotic  (/  =  0.95  Hz) 
trajectory.  These  two  frequencies  correspond  to  Fig.  7 
and  Fig.  9  respectively.  This  highlights  the  drastic  dif¬ 
ference  between  the  peak  counts  for  periodic  and  chaotic 
trajectories.  The  peak  counts  fall  off  very  rapidly  with 
threshold  for  the  periodic  case,  while  they  only  gradu¬ 
ally  drop  off  for  the  chaotic  case.  This  shows  that  the 
output  (binary  output  of  either  chaotic  or  nonchaotic)  is 
quite  insensitive  to  the  selection  of  the  cut-off  threshold. 
The  method  is  still  likely  limited  to  systems  with  rela¬ 
tively  high  signal  to  noise  ratios  where  the  noise  does  not 
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FIG.  10.  Spectrograms,  (a)  numerical,  (b)  experimental. 


swamp  the  frequency  spectra  of  the  response.  Of  course, 
the  presence  of  noise  also  makes  the  determination  of 
LE’s  more  difficult.  However,  the  system  noise  level  is 
not  atypical  of  periodically  forced  low-order  mechanical 
systems. 

In  order  to  verify  the  above  bifurcation  diagrams,  sev¬ 
eral  other  methods  were  also  applied.  Figure  13  shows 
numerical  and  experimental  one-parameter  amplitude  bi¬ 
furcation  diagrams  over  the  same  frequency  range.  The 
dashed  vertical  lines  again  show  the  locations  of  the  four 
samples  in  Fig.’s  6  to  9.  The  ordinate  axis  shows  the 
amplitude  difference  between  a  peak  and  the  following 
valley  in  the  time  series  of  0(t).  To  develop  this  plot 
an  individual  simulation/test  was  run  for  each  forcing 
frequency  (400  simulations  numerically,  58  tests  experi¬ 
mentally).  The  system  was  again  allowed  to  run  for  an 
initial  period  of  time  to  allow  transients  to  decay.  For  a 
periodic  time  series  all  of  the  data  points  will  lay  on  top 


FIG.  11.  Plot  of  peak  counts  of  DFT,  threshold  =  1%.  Nu¬ 
merical  (red),  experimental  (black). 
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FIG.  12.  Sensitivity  of  peak  counts  to  cut-off  threshold. 
Threshold  =  0%  (black),  1%  (blue),  5%  (red),  10%  (green). 


FIG.  13.  Amplitude  bifurcation  diagram,  (a)  numerical,  (b) 
experimental. 


Figure  14  shows  the  numerically  obtained  largest  LE 
bifurcation  diagram  for  the  same  parameter  range.  The 
curve  was  obtained  by  simulating  the  equation  of  mo¬ 
tion,  and  calculating  the  largest  LE  using  the  method 
in5 .  Once  again  the  dashed  vertical  lines  denote  the  four 
samples  in  Fig.’s  6  to  9.  The  regions  of  positive  largest 
LE  agree  well  with  the  regions  of  chaos  in  Fig.’s  10  to  13. 
Note  that  for  a  nonchaotic  trajectory  the  true  largest  LE 
of  a  temporally  forced  system  is  zero  since  the  time  state 
is  non-converging.  However  since  the  numerical  method 
has  direct  access  to  each  state  one  may  remove  time  from 
the  divergence  calculations  and  therefore  obtain  negative 
LE’s. 


of  one  another,  while  for  a  chaotic  time  series  they  will 
be  scattered  about.  It  is  interesting  in  the  chaotic  band 
between  approximately  /  =  0.8 Hz  and  /  =  1.0 Hz  that 
there  is  a  strip  of  mid-range  Ad  values  that  are  never 
visited.  This  makes  sense  since  in  a  system  with  a  pe¬ 
riodic  potential  the  oscillations  will  either  be  within  a 
single  potential  well  and  therefore  small,  or  crossing  be¬ 
tween  the  two  wells  and  therefore  large.  The  agreement 
between  the  experimental  and  numerical  results  is  excel¬ 
lent  and  they  both  show  a  resonant  jump  in  amplitude 
near  /  =  0.65 Hz  which  is  typical  of  ’softening’  spring 
systems23.  There  is  also  excellent  agreement  with  the 
peak  count  plots. 


FIG.  14.  Numerical  LE  bifurcation  diagram. 

Determination  of  LE’s  from  experimental  data  is  a  sen¬ 
sitive  procedure.  In  calculating  LE’s  numerically  from  an 
equation  of  motion  one  has  the  ability  to  either  normal¬ 
ize  and  restart  the  analysis  of  divergence  rates  of  nearby 
trajectories  or  to  simply  linearize  the  divergence  rates. 
This  makes  it  possible  to  calculate  the  LE  directly  by 
its  definition  as  an  exponential  divergence  rate.  This  is 
not  possible  with  experimental  data.  One  useful  method 
however,  especially  with  small  data  sets,  is  discussed  in' . 
In  this  method  the  phase  portrait  is  first  reconstructed 
and  then  the  nearest  neighbor  is  found  for  each  data 
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point.  Then  the  pairs  of  nearest  neighbors  are  stepped 
forward  in  time  and  their  divergence  rate  is  monitored. 
Finally  the  average  of  the  logarithm  of  the  divergence 
rate  of  nearest  neighbors  is  taken  over  all  sets  of  pairs 
at  each  point  in  time  (relative  time  for  each  pair).  In 
the  best  case  there  is  a  relatively  long  lasting  period  of 
time  in  which  the  plot  of  the  logarithm  of  the  divergence 
rate  is  approximately  linear.  The  slope  of  the  linear  por¬ 
tion  is  of  course  the  divergence  rate,  and  therefore  the 
LE  of  the  system.  Of  course,  as  two  nearby  trajectories 
grow  far  apart  the  dynamics  of  their  divergence  rate  will 
become  more  and  more  nonlinear  and  hence  the  loga¬ 
rithm  of  the  divergence  will  eventually  lose  its  linearity. 
The  time  scale  over  which  the  linear  behaviour  occurs 
and  the  largest  LE  may  be  calculated  is  highly  system 
and  parameter  dependent  making  it  very  difficult  to  tune 
for  a  hands-off  analysis  over  parameter  space.  For  many 
systems,  especially  autonomous  or  low  dimensional  sys¬ 
tems,  the  nearest  neighbors  divergence  approach  works 
very  well  for  calculating  LE’s  (See  Appendix).  However 
for  the  forced  oscillator  studied  herein  the  results  were 
less  impressive.  Figure  15  shows  the  logarithm  of  the  di¬ 
vergence  rate  of  initially  nearby  trajectories  as  a  function 
of  time  for  two  different  sets  of  forcing  parameters.  Part 
(a)  is  for  the  chaotic  response  in  Fig.  9  while  part  (b)  is 
for  the  periodic  response  in  7.  The  red  lines  (left  vertical 
axes)  shows  the  average  of  the  Logarithm  of  the  diver¬ 
gence  of  nearest  neighbors  for  experimental  data,  while 
the  blue  lines  (right  vertical  axes)  shows  the  same  (using 
the  same  method)  on  numerically  simulated  data  from 
the  equation  of  motion.  The  experimental  curves  have 
a  much  larger  y-intercept  as  the  nearest  neighbors  are 
much  farther  apart  with  experimental  data.  The  slope 
of  the  dashed  lines  is  equal  to  the  exact  LE’s,  that  is  the 
LE  obtained  directly  from  the  equation  of  motion  earlier. 
This  dashed  line  was  placed  artificially  and  its  y-intercept 
is  not  particularly  meaningful.  The  time  lag  (for  this  case 
it  was  1  /4  of  the  forcing  period)  which  is  used  to  recon¬ 
struct  the  phase  portrait  can  effect  the  result.  However 
the  curves  in  Fig.  15  were  relatively  insensitive  to  this 
parameter,  a  fact  which  is  also  mentioned  in7.  It  is  clear 
from  part  (a)  that  tuning  the  time  scale  over  which  to 
calculate  the  LE  would  be  difficult  for  this  system  as 
there  is  no  clear  linear  portion  of  the  line.  For  part(b) 
one  would  hope  that  the  experimentally  obtained  LE  for 
a  periodic  trajectory  would  be  either  negative,  zero,  or 
at  least  much  smaller  than  those  for  a  chaotic  response. 
This  is  the  case  for  the  blue  line,  but  the  experimental 
data  still  appears  to  have  a  positive  slope.  It  is  likely 
that  this  is  due  to  the  real  noise  in  the  system,  as  the 
phase  portrait  in  Fig.  7  is  clearly  periodic. 

In  general  it  is  not  immediately  clear  whether  the 
heuristic  method  or  the  LE  approach  require  longer  time 
series  to  obtain  a  good  result,  however  for  this  system  it 
appears  that  the  heuristic  method  yields  better  results. 
The  experimental  time  series  used  in  this  system  were 
all  100s  long,  which  is  approximately  16  forcing  cycles 
in  the  vicinity  of  the  1  Hz  forcing  region  of  the  bifurca¬ 


tion  diagrams.  This  was  done  because  a  large  number  of 
parameters  were  investigated.  To  the  naked  eye  it  is  rel¬ 
atively  easy  to  distinguish  between  periodic  and  chaotic 
responses  over  this  number  of  cycles.  This  is  also  usually 
enough  cycles  to  capture  quasi-periodicity  on  the  DFT 
of  the  response,  as  long  the  Nyquist  frequency  is  high 
enough  to  capture  all  of  the  dominant  frequencies  (the 
sampling  rate  used  was  200  Hz  which  was  well  above 
any  contributing  frequencies  in  the  response).  However, 
especially  for  a  chaotic  response,  this  may  not  provide 
enough  pairs  of  nearby  neighbors  with  which  to  yield  a 
good  approximation  of  the  LE. 


FIG.  15.  Average  of  logarithm  of  divergence  of  nearest  neigh¬ 
bors,  L  =  5.5 cm,  (a)  /  =  0.95 Hz,  (b)  /  =  0.55 Hz. 

In  order  to  obtain  the  boundaries  between  chaotic  and 
periodic  response  in  two-parameter  space,  (/,  L ),  one 
can  plot  the  peak  count  as  a  contour  plot  versus  any  two 
parameters.  Figure  16  shows  the  results  from  numerical 
LE  calculations  where  regions  colored  orange  indicate  pa¬ 
rameters  with  a  chaotic  response  ( LE  >  0)  while  regions 
in  red  indicate  a  periodic  response  ( LE  <  0).  The  hor¬ 
izontal  line  at  L  =  5.5 cm  denotes  the  forcing  amplitude 
of  the  one-parameter  bifurcation  diagrams  presented  ear¬ 
lier.  The  resolution  of  the  plot  is  200  points  in  the  fre¬ 
quency  axis  and  200  in  the  forcing  amplitude  axis  yielding 
40,000  total  points. 

The  two-parameter  plot  is  where  the  heuristic  method 
is  especially  useful  over  the  standard  spectrogram.  Fig¬ 
ure  17  shows  the  results  of  the  two-parameter  investi¬ 
gation  using  the  heuristic  method  on  both  numerically 
simulated  and  experimental  data.  The  peak  count  cut-off 
threshold  used  was  1%  and  5%  respectively  for  the  nu¬ 
merical  and  experimental  cases.  Selection  of  these  values 
was  aided  by  referencing  Fig.’s  11  and  12.  The  numerical 
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threshold  could  be  chosen  smaller  since  the  only  source 
of  ’noise’  is  due  to  the  numerical  errors.  To  allow  for  a  di¬ 
rect  comparison  with  LE’s  it  was  decided  to  set  a  second 
threshold  peak  count  number  separating  classification  as 
either  chaotic  or  periodic.  This  was  made  relatively  easy 
due  to  the  clear  dichotomy  in  the  peak  counts  between 
chaotic  and  periodic  responses  shown  earlier.  This  second 
threshold  was  chosen  to  be  40  and  20  peaks  respectively 
for  the  numerical  and  experimental  cases.  These  values 
were  chosen  as  points  approximately  mid-way  between 
typical  chaotic  and  periodic  response  peak  counts  by  ref¬ 
erencing  the  one-parameter  peak  count  plots  in  Fig.’s  14 
and  15.  The  resolution  for  part  (a)  is  40,000  data  points 
(same  as  Fig.  16)  while  the  resolution  for  the  experimen¬ 
tally  obtained  plot  in  part  b)  is  only  58x15.  The  agree¬ 
ment  is  once  again  excellent,  as  is  the  agreement  between 
the  LE’s  and  the  heuristic  method.  It  is  important  to 
note  that  the  standard  LE  approach  also  frequently  re¬ 
quires  calibration  via  visual  inspection  for  experimental 
data.  The  run  time  over  which  to  determine  the  LE  is 
especially  difficult  to  perform  automatically  and  is  fre¬ 
quently  somewhat  ambiguous  as  seen  in  Fig.  15. 
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FIG.  17.  Chaotic  boundaries  in  parameter  space:  heuristic 
method,  (a)  Numerical,  (b)  experimental. 
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IV.  CONCLUSIONS 


A  heuristic  method  for  identifying  chaos  is  developed 
and  tested  both  numerically  and  experimentally  on  a 
nonlinear  mechanical  oscillator.  The  method  is  used  to 
create  several  bifurcation  diagrams  and  shows  excellent 
agreement  with  numerically  obtained  Lyapunov  expo¬ 
nents,  much  better  in  fact  than  the  experimentally  ob¬ 
tained  Lyapunov  exponents. 
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Appendix:  The  Lorenz  Equations 

In  order  to  show  the  applicability  of  the  method  to 
other  systems  it  has  been  applied  numerically  to  the  fa¬ 
miliar  Lorenz  equations24,25 

x  =  a(y-  x), 

y  =  rx  —  y  —  xz,  (A.l) 

z  =  xy  —  bz, 

with  parameters  (it,  r,  6)  and  states  (x,y,z).  Figure  18 
shows  a  typical  periodic  and  chaotic  response  for  the 
Lorenz  equations.  Figure  19  shows  the  chaotic  bound¬ 
aries  in  (cr,  r)  space  with  b  =  8/3  obtained  using  the 
heuristic  method  and  the  using  LE’s.  The  white  dots 
shown  in  the  figure  denote  the  location  of  the  samples 
shown  in  Fig.  18.  The  heuristic  method  was  carried  out 
using  a  peak  count  cut-off  threshold  of  1%  of  the  maxi¬ 
mum  peak. 
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FIG.  18.  Sample  responses  (numerical)  from  the  Lorenz  equa¬ 
tions,  r=200,  b=8/3,  (a)  a  =  100,  (b)  a  =  25. 
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FIG.  19.  Chaotic  boundaries  in  two-parameter  space  for  the 
Lorenz  system  for  6  =  8/3.  (a)  Heuristic  method,  (b)  LE’s. 


The  experimental  nearest  neighbors  divergence 
method  works  better  on  the  Lorenz  system  than  it  does 
on  the  forced  oscillator  (Fig.  12).  Figure  20  shows  the 
results  of  the  nearest  neighbors  approach  versus  direct 
calculation  of  LE’s.  These  are  ID  slices  as  indicated 
by  the  horizontal  dashed  lines  in  19,  i.e.  three  chaotic 
zones  when  r  =  500,  and  two  when  r  =  200.  Once 
again  the  nearest  neighbors  method  here  is  applied 
to  data  simulated  from  the  equation  of  motion,  and 
then  compared  to  the  direct  method  of  LE  calculation. 
Parts  (a)  and  (b)  show  the  average  logarithm  of  the 
divergence  of  nearest  neighbors  for  two  different  sets  of 
parameters  located  in  two  different  bands  of  chaos  for 
r  =  500.  The  time  scale  over  which  the  LE  may  be 
calculated  (the  linear  portion)  is  much  easier  to  identify 


and  appears  to  be  much  more  stable  for  the  Lorenz 
system.  The  dashed  lines  superimposed  on  each  plot  has 
the  slope  of  the  LE  which  was  obtained  by  the  direct 
numerical  approach.  Part  (c)  shows  a  comparison  of 
the  one-parameter  LE  bifurcation  plots  obtained  using 
the  experimental  approach  (black  line)  and  the  direct 
numerical  approach  (red  line)  for  r  =  500,  where  the 
two  vertical  dotted  lines  indicate  the  location  of  parts 
(a)  and  (b).  Finally  part  (d)  shows  the  same  comparison 
of  the  bifurcation  diagrams,  but  this  time  for  r  =  200, 
while  the  two  vertical  dashed  lines  now  denote  the 
location  of  the  phase  portraits  in  Fig.  18.  For  the 
experimental  data  approach  the  LE  was  calculated  by 
obtaining  the  approximate  slope  of  the  logarithm  of 
the  divergence  plot  from  t  =  0s  to  t  =  0.3s  as  this 
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appears  to  be  reasonably  a  stable  zone  where  the  line 
remains  nearly  linear.  The  curves  were  all  obtained 
with  50s  simulations  and  5000  data  points  and  the  time 
lag  used  for  phase  space  reconstruction  was  11  time 
steps.  These  are  all  the  same  parameters  that  were  used 
in7.  The  difference  in  the  magnitude  between  the  two 
methods  could  likely  be  reduced  by  a  more  thorough 
calibration  of  the  time  lag  and  slope  measurement  time 
span,  however  the  agreement  on  the  locations  of  the 
bifurcations  between  periodic  and  chaotic  behaviours  is 
excellent. 


i  i 


FIG.  20.  Comparison  of  LE  calculation  methods  for  b  = 
8/3, r  =  500  for  (a),(b),(c),  and  r  =  200  for  (d).  Parts  (a) 
a  =  50,  and  (b)  a  =  130,  show  the  average  logarithm  of  the 
divergence  rate  of  nearest  neighbors.  Parts  (c)  and  (d)  show 
bifurcation  diagrams  from  the  direct  numerical  method  (red) , 
and  the  experimental  approach  (black). 
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A  non-linear  finite  element  formulation  (three  dimensional  continuum  elements)  is  implemented  and 
used  for  modeling  dynamic  snap-through  in  beams  with  initial  curvature.  We  identify  a  non-trivial 
(non-flat)  configuration  of  the  beam  at  a  critical  temperature  value  below  which  the  beam  will  no 
longer  experience  snap-through  under  any  magnitude  of  applied  quasi-static  load  for  beams  with 
various  curvatures.  The  critical  temperature  is  shown  to  successfully  eliminate  snap-through  in  dynamic 
simulations  at  quasistatic  loading  rates.  Thermomechanical  coupling  is  included  in  order  to  model  a 
physically  minimal  amount  of  damping  in  the  system,  and  the  resulting  post-snap  vibrations  are  shown 
to  be  thermoelastically  damped.  We  propose  a  test  to  determine  the  critical  snap-free  temperature  for 
members  of  general  geometry  and  loading  pattern;  the  analogy  between  mechanical  prestress  and 
thermal  strain  that  holds  between  the  static  and  dynamic  simulations  is  used  to  suggest  a  simple 
method  for  reducing  the  vulnerability  of  thin-walled  structural  members  to  dynamic  snap-through  in 
members  of  large  initial  curvature  via  the  introduction  of  initial  pretension. 
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1.  Introduction 

Thin  structural  members  can  dynamically  jump  between 
multiple  equilibrium  configurations  when  subjected  to  mechan¬ 
ical  forces,  acoustic  vibrations,  and  thermal  loads  like  those 
encountered  by  aerospace  vehicles  in  extreme  operating  condi¬ 
tions.  This  process,  commonly  referred  to  as  snap-through,  can 
cause  large-amplitude  structural  vibrations,  induce  fatigue,  and 
lead  to  global  instability.  Such  vibrations  can  have  a  chaotic 
pattern  and  therefore  can  be  difficult  to  control  [1].  Avoiding 
snap-through  is  therefore  highly  desirable  in  the  design  of  thin- 
walled  structures  such  as  aircraft  and  spacecraft. 

In  this  paper,  the  existence  of  a  critical  temperature  below 
which  an  initially  curved  beam  will  no  longer  experience  snap- 
through  at  any  applied  quasi-static  load  level  is  demonstrated. 
This  temperature  provides  a  lower  bound  below  which  snap- 
through  instability  no  longer  occurs.  This  temperature  also 
corresponds  to  a  non-trivial  (that  is,  non-flat)  deformed  beam 
configuration,  a  fact  that  is  not  apparent  through  methods  used  to 
characterize  snap-through  that  are  limited  to  modeling  small 
deflections.  We  argue  that  such  a  limit  can  be  obtained  for  a 
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variety  of  static  loads  and  that  it  also  reduces  the  amplitude  of  the 
oscillations  when  loads  are  applied  at  dynamic  rates.  Through 
numerical  experiments  we  show  this  limit  to  be  found  at  states 
characterized  by  specific  properties  associated  with  the  energy  of 
the  system  and  the  directional  derivatives  of  this  energy  function. 

Experimental  characterizations  of  snap-through  under  mechan¬ 
ical  vibrations  and  combined  thermal  and  mechanical  loads  have 
existed  for  several  decades  [2-4],  but  numerical  modeling  of 
snap-through  remains  an  area  of  active  research.  In  the  general 
case,  snap-through  involves  combined  flexural,  shear,  and  normal 
stresses,  large  time  varying  deflections  and  rotations,  and  thermal 
effects.  Non-linear  coupling  between  these  effects  makes  snap- 
through  challenging  to  model  accurately.  An  extensive  literature 
exists  on  the  subject,  with  models  describing  different  degrees  of 
geometric  non-linearity  [5,6],  material  non-linearity  [7-9],  thermal 
and  acoustic  loading  [10-12],  and  various  combinations  of  these 
factors  [13,14], 

Previous  research  on  snap-through  falls  into  three  general 
categories:  (1)  limited-deflection  models  (based  on  the  von 
Karmann  or  Duffing  equations  for  plates),  (2)  elastica  models 
(based  on  specialized  analytical  techniques  going  back  to  Euler), 
and  (3)  non-linear  (finite-deformation)  finite  element  models. 

The  first  category  suffers  from  the  most  significant  limitation: 
these  models  can  only  describe  member  deflections  of  less  than 
2.5  times  the  member  thickness  [5[.  This  is  especially  problematic 
because  large  initial  curvatures  that  may  be  present  by  design, 
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e.g.,  in  a  rocket  booster  hull  or  aircraft  wing,  will  lead  to  larger 
snap-through  loads  and  therefore  more  severe  and  damaging 
post-snap  vibrations,  none  of  which  can  be  modeled  with  limited 
deflection  formulations. 

Elastica  methods  are  valid  at  a  theoretically  unlimited  range  of 
deflections  and  initial  curvatures,  at  the  cost  of  increased  math¬ 
ematical  complexity.  These  techniques  have  only  recently  been 
extended  to  the  solution  of  thermal  and  snap-through  problems 
[15-17].  This  is  most  likely  due  to  the  increased  ease  of  solving 
the  equations  numerically,  e.g.,  by  shooting  methods  [18,19]. 
Although  there  is  no  apparent  limitation  to  obtaining  the  results 
discussed  in  this  paper  by  elastica  methods  (other  than  possible 
limitations  to  particular  instances  of  the  theory,  e.g.,  the  common 
omission  of  shear  strain),  they  are  not  as  well-developed  for 
general  classes  of  problems  as  non-linear  finite  element  methods. 

Non-linear  finite  element  methods  are  the  most  general  of  the 
methods  available  to  model  snap-through,  an  advantage  for  which 
they  sacrifice  the  direct,  problem-specific  analytical  insight  that  is 
available  via  the  first  two  methods.  This  lack  of  transparency  is 
more  than  compensated  by  the  versatility  to  model  geometric 
non-linearity,  thermomechanical  coupling,  large  strains,  material 
non-linearity,  fluid-structure  interaction  in  aircraft,  etc. 

The  non-linear  finite  element  method  is  chosen  for  this  study,  it  is 
capable  of  including  the  full  range  of  non-linearities  that  influ¬ 
ence  snap-through  and  can  easily  model  large  deflections.  Beam 
and  shell  models,  usually  the  most  computationally  efficient 
choice  for  finite  element  modeling  of  thin  structural  members, 
always  involve  some  simplifications  of  the  underlying  three- 
dimensional  kinematics  that  can  lead  to  artificial  stiffness  under 
particular  load  states  (locking).  These  kinematic  assumptions  that 
are  built  into  the  formulation  of  structural  elements  can  lead  to 
inaccurate  solutions  [20],  In  this  work  we  avoid  such  issues  as 
well  as  the  locking  sometimes  present  when  linear  elements  are 
used  by  choosing  to  work  with  three-dimensional  (solid)  quad¬ 
ratic  elements  only.  In  this  paper,  the  numerical  simulations  are 
performed  with  the  Finite  Element  Analysis  Program  (FEAP),  an 
open  source  research  code,  which  provides  a  framework  for  finite 
element  simulations  where  we  can  formulate  and  implement 
additional  elements,  constitutive  models  and  solution  schemes 
via  user  subroutines.  The  formulations  utilized  for  the  analysis 
discussed  in  this  paper  are  a  mix  of  FEAP  original  elements  and 
user  routines  [21]. 

By  thermomechanical  coupling  we  refer  only  to  the  coupling 
between  elastic  deformation  and  thermal  effects  via  the  thermal 
strain  terms  in  the  equations  of  mechanical  equilibrium  and  the 
structural  elastic  heating  term  in  the  heat  equation.  Although 
thermoelastic  coupling  is  present  in  all  materials,  it  is  often 
neglected;  however,  for  large-amplitude  vibrations,  the  thermal 
gradient  induced  between  the  compressive  and  tensile  fibers  of  a 
vibrating  member  due  to  purely  elastic  deformation  results  in 
significant  heat  conduction  and  therefore  loss  of  energy  via 
thermoelastic  damping.  We  do  not  include  any  other  mechanisms 
of  dissipation  besides  thermoelastic  damping,  so  the  damping  in 
our  study  is  physically  minimal,  and  the  resulting  vibrations  are 
exaggerated  relative  to  real  physical  systems. 

In  many  studies  [22,23],  thermal  strain  due  to  applied  tem¬ 
perature  is  included  without  modeling  thermoelastic  coupling: 
temperature  changes  are  thus  treated  as  a  simple  mechanical 
expansion  or  contraction  of  the  material.  This  can  be  useful  in 
tracing  static  solution  paths,  and  may  also  give  information  about 
the  temperature-sensitivity  of  a  particular  structure,  but  the 
absence  of  thermoelastic  coupling  can  impact  the  accuracy  of 
dynamic  simulations,  especially  as  amplitudes  of  vibration 
become  large.  We  employ  the  strain-only  method  in  our  static 
simulations,  but  move  to  full  coupling  in  order  to  capture  relevant 
dynamic  coupling  effects. 


Modeling  of  thermomechanical  coupling  also  requires  the 
selection  of  an  appropriate  thermomechanical  material  constitu¬ 
tive  law.  A  hyperelastic  constitutive  law  based  on  a  modified 
neo-Flookean  material  that  is  appropriate  for  metals  has  been 
proposed  in  previous  studies  [24].  Since  we  will  be  primarily 
concerned  with  metals,  the  modified  thermomechanical  neo- 
Hookean  material  law  is  appropriate  and  will  be  valid  at  a  much 
larger  range  of  strains  than  the  thermomechanical  St.  Venant- 
Kirchhoff  law. 

Failure  to  account  for  any  of  the  above  modeling  concerns  can 
significantly  impact  the  accuracy  of  dynamic  simulation  results. 
The  formulation  adopted  for  this  study  accounts  for  the  geometric 
non-linearity,  large  strains,  and  thermomechanical  coupling 
effects  required  to  describe  snap-through.  It  is  a  non-linear  finite 
element  formulation  that  uses  the  adiabatic  staggered  scheme 
and  the  constitutive  law  detailed  in  [24]  and  quadratic  hexahe- 
dral  (or  tetrahedral)  elements. 

The  initial  boundary  value  problem  is  formulated  as  follows: 
For  all  tsl  find  the  motion  (</>)  and  temperature  (T)  fields 
such  that 

a2 

Poj»2<l>  =  divP+b, 

ot 

cT  =  V-K.-f  div[q/J]+Tl,  (1) 

with  boundary  conditions  Pn0  =  t  on  f„  x  I,  </>  =  </>  on  F^,  x  D, 
qn0  =  q  on  Fq  x  D,  T=  T  on  Fq  x  1  and  initial  conditions  <j)\t  =  0  =  I 
in  Q,  (d/8t)tf>\t  =  0  =  V0  in  Q,  T|t  =  0  =  T0  in  Q.  12c  R"11  is  the 
domain  in  the  reference  configuration,  nd  is  the  number  of  spatial 
dimensions  in  the  problem,  p0  is  the  reference  configuration 
density,  P  is  the  first  Piola-Kirchhoff  stress  tensor,  b  is  the 
prescribed  body  force,  c  is  the  heat  capacity,  D  is  the  mechanical 
dissipation,  K.  is  the  heating  from  the  Joule  effect  [25],  Tl  is  the 
prescribed  heat  source  term,  and  J  is  the  Jacobian  of  the  deforma¬ 
tion.  •  denotes  the  prescribed  value  of  the  quantity  •  over  the 
appropriate  boundary  region  or  at  the  initial  time.  The  generic 
problem  described  here  can  be  completed  with  a  constitutive  law 
prevailing  in  the  body  and  the  Fourier  law  for  heat  conduction  is 
assumed  to  relate  the  local  heat  flux  q  to  the  temperature 
gradient,  q=-[/d]VT,  where  k  is  the  thermal  conductivity  and 
1  is  the  identity  tensor.  For  the  purpose  of  this  paper,  the 
constitutive  law  is  assumed  to  have  temperature-dependent 
constitutive  moduli.  After  a  spatial  discretization  is  applied  (e.g., 
finite  element)  the  system  (1)  can  be  expressed  as  a  system  of 
ordinary  differential  equations  with  two  coupled  partitions: 
mechanical  (second  order)  and  thermal  (first  order). 

For  the  transient  solution  we  use  the  adiabatic  staggered 
scheme  for  coupled  thermoelastic  boundary  value  problems 
developed  in  [24];  this  approach  splits  the  solution  into  partitions 
such  that  the  dissipative  property  of  the  original  problem  is 
maintained  in  the  partitioned  problem.  The  scheme  consists  of 
(1)  a  standard  finite  deformation  mechanical  phase,  formulated 
such  that  it  is  solved  at  constant  entropy,  and  (2)  a  heat  conduc¬ 
tion  phase,  formulated  such  that  it  is  solved  at  fixed  deformation. 

We  also  adopted  the  constitutive  law  from  [24,  p.  760, 
Eqs.  (82)-(83)[,  for  a  regularized  compressible  neo-Hookean 
material.  This  is  a  law  that  performs  well  for  metals  and  is  given 
by  the  thermomechanical  strain  energy  density  function: 

<F(C,0)  =  W(C)  +  U(J)+M(J,0)  +  T(0),  (2) 

where  W,  U,  M,  and  T  are  the  mechanical  deviatoric,  mechanical 
volumetric,  thermomechanical  coupling,  and  thermal-only  terms 
respectively,  given  by 

W(C)=i  /i(tr(C)-3)  =  i/i(r2/3  tr(C)-3), 

U(J)  =  lK(lnJ)2, 
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M(J,& )  =  -3Koc(0-0ref)ln  J, 


T(0)  =  pcm 


( 0-0ret)-0  In 


where  C  is  the  right  Cauchy-Green  tensor,  /( is  the  shear  modulus, 
cm  is  the  mass  heat  capacity,  0ref  is  the  reference  temperature, 
I<  is  the  bulk  modulus,  a  is  the  coefficient  of  thermal  expansion 
and  J  is  the  Jacobian  of  deformation. 

For  the  detailed  derivation  of  the  components  of  the  finite 
element  formulations  utilized  in  this  work,  the  interested  reader 
is  referred  to  [26]. 

The  rest  of  the  paper  is  organized  as  follows.  In  Section  2  we 
describe  the  geometry  of  the  system  under  consideration  and 
briefly  describe  the  solution  methods  utilized  to  retrieve  equili¬ 
brium  paths.  Section  3  is  dedicated  to  numerical  experiments  that 
trace  such  equilibrium  paths.  In  this  section,  we  identify  special 
equilibrium  configurations  corresponding  to  boundaries  in  the 
parameter  space  that  separate  domains  with  different  stability 
behaviors.  Among  these,  the  quasi-static  temperature  limit  that 
we  will  refer  to  as  monotonic,  i.e.,  that  value  of  the  temperature 
below  which  the  equilibrium  path  is  a  monotonic  curve,  and 
snap-through  is  no  longer  encountered.  Section  4  demonstrates, 
through  a  variety  of  transient  simulations,  that  the  monotonic 
temperature  also  identifies  a  limit  below  which  the  amplitude  of 
the  dynamic  oscillations  is  drastically  reduced  and  asymptotically 
reaches  0  at  temperatures  in  the  neighborhood  of  the  quasi-static 
limit.  Section  5  provides  a  generalization  of  the  concept  of 
monotonic  temperature  and  shows  that  information  regarding 
such  stability  limits  is  no  longer  available  when  some  specific 
finite  element  formulations  (e.g.,  the  Timoshenko  beam)  are 
utilized.  When  kinematic  approximations  such  as  those  used  by 
the  Timoshenko  beam  are  used,  information  available  in  the 
higher  derivatives  of  the  energy  is  lost  in  such  formulations. 


2.  Static  analysis  of  a  curved  beam 

A  planar  curved  beam  described  by  an  arc  of  a  circle  is  chosen 
as  the  test  problem  to  explore  the  effect  of  applied  temperature 
and  initial  beam  curvature  on  snap-through.  The  problem  is 
solved  in  two  phases  (mechanical  and  thermal).  The  beam 
geometry  is  completely  specified  by  the  length  of  the  horizontal 
projection  between  the  end  points  L,  the  projection  of  the  lengths 
of  the  supports  LBC,  and  the  radius  of  curvature  of  the  beam  R 
(Fig.  1).  The  primary  load  case  considered  is  a  point  load  P  in  the 
negative  y-direction  located  at  the  midpoint  between  the  sup¬ 
ports  (Fig.  2a).  A  distributed  load  p,  applied  as  a  uniform  force  in 
the  negative  y-direction  on  each  node  of  the  finite  element  mesh, 
is  also  considered  (Fig.  2b). 

In  simulations,  eight  different  beams  were  utilized,  with  R 
varying  from  762  mm  (30  in.)  to  5080  mm  (200  in.).  The  effect  of 
temperature  variation  on  the  load-deflection  behavior  in  the  case 
of  larger  curvatures  was  found  to  be  very  small,  so  such  test 
problems  were  abandoned. 

Table  1  summarizes  the  geometry  of  the  beams:  M  is  the  arch 
rise,  k  =  1  /R  is  the  curvature  of  the  beam,  9  is  the  angle  subtended 
by  the  beam,  and  M/h  is  the  ratio  of  M  to  the  thickness  h  of  the 
beam.  Fig.  3  shows  the  different  initial  beam  curvatures.  For  all 
beams,  the  projection  length  is  L=304.8  mm  (12  in.),  the  thick¬ 
ness  li  =  0.508  mm  (0.02  in.)  and  the  transverse  depth  b=  12.7  mm 
(0.5  in.).  The  ends  of  the  beam  are  held  fixed  over  a  length 
Ibc=28.7  mm  (1.125  in.)  on  each  end.  Material  properties  are 
those  of  steel,  given  in  Table  2. 

In  this  paper  we  concentrate  our  attention  to  a  system 
that  is  effectively  two-dimensional;  nevertheless,  we  use  a  finite 
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Fig.  2.  Concentrated  (a)  and  distributed  (b)  loads.  The  point  load  is  the  primary 
case  considered. 


Table  1 

Geometry  of  curved  beams. 


Beam 

R  (mm) 

M  (mm) 

k  (1/mm) 

6/2  (rad) 

M/h 

1 

762.0 

21.84 

1.312  x  1(T3 

0.2400 

86.0 

2 

1270.0 

12.99 

0.787  xKT3 

0.1431 

51.1 

3 

1828.8 

8.98 

0.546  x  10~3 

0.0990 

35.3 

4 

2133.6 

7.70 

0.469  x  10~3 

0.0850 

30.3 

5 

2438.4 

6.74 

0.410  x  10-3 

0.0744 

26.5 

6 

3048.0 

5.39 

0.328  xlO-3 

0.0595 

21.2 

7 

3810.0 

4.31 

0.262  x  10~3 

0.0475 

17.0 

8 

5080.0 

3.23 

0.169  x  10~3 

0.0357 

12.7 

element  model  with  full  three-dimensional  capabilities.  Mesh 
refinement  studies  were  performed  and  the  level  of  refinement 
adequate  to  accurately  capture  the  snap-through  loads  lead 
to  a  model  with  approximately  20,000  degrees  of  freedom.  Note 
that  linear  elements  suffer  from  shear  locking  and  are  not 
appropriate  for  snap-through  problems,  or  for  any  problems 
involving  large  deformation.  This  problem  is  avoided  by  using 
quadratic  elements. 

The  formulation  adopted  [24]  treats  the  coupling  through  an 
adiabatic  staggered  approach.  For  the  quasistatic  simulations 
of  systems  under  constant  temperature  presented  in  the  next 
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Fig.  3.  Initial  beam  configurations. 


Table  2 

Material  properties  of  thermoelastic  beam. 


Property 

Symbol 

Value 

Units 

Young’s  modulus 

E 

206,483 

N/mm2 

Poisson’s  ratio 

V 

0.28 

- 

Density 

p 

7.834  x  10-9 

N  s2/mm4 

Conductivity 

k 

45.0 

N/sI< 

Specific  heat 

Cm 

to 

O 

X 

ro 

mm2/s2  K 

Thermal  expansion 

a 

14  x 10~6 

mm/mm  K 

section,  the  mechanical  phase  is  solved  in  isolation  from  the 
thermal  phase,  treating  the  applied  temperature  change  as  a 
conservative,  purely  mechanical  expansion  or  contraction  of  the 
beam  material.  Many  other  studies  have  employed  this  method, 
but,  with  a  few  exceptions  [15,16,6],  they  have  been  limited  to 
beams  of  small  initial  curvature. 

The  case  of  the  point  load  along  the  line  of  symmetry  of  the 
beam  provides  an  easy  interpretation  of  the  simulation  results. 
If  the  symmetry  is  not  broken  by  a  bifurcation,  as  may  occur  if  the 
second,  asymmetric  buckling  mode  emerges  [27],  the  plots  of 
midpoint  deflection  d  versus  load  value  P  will  indicate  loss  of 
stability  and  locate  the  snap-through  and  snap-back  loads  at  the 
points  where  SP/dd  =  0.  The  ability  to  obtain  such  information  is 
completely  dependent  on  the  symmetry  of  the  problem:  in  cases 
without  a  clear  line  of  symmetry  it  is  no  longer  possible  to  obtain 
information  about  the  snap-through  point  from  the  load-deflec¬ 
tion  curve,  and  even  in  the  symmetric  case  of  the  distributed  load 
this  breaks  down  at  heightened  temperatures.  The  simplicity  of 
the  symmetric  concentrated  load  case  will  make  it  easier  to 
develop  insight  into  the  stability  behavior  of  the  beam,  which 
we  will  then  need  to  generalize  to  cases  where  important 
behaviors  are  no  longer  apparent  from  the  load-deflection  curves. 

2.2.  Static  solution  paths 

Snap-through  is  inherently  a  dynamic  phenomenon;  since  the 
static  equilibrium  path  is  interrupted  by  an  unstable  region,  as 
the  load  is  increased,  the  system  must  dynamically  jump  past  the 
unstable  region  and  onto  a  stable  region  capable  of  bearing  loads 
above  the  snap-through  load.  Unlike  column  buckling,  there  is  no 
stable  branch  that  the  system  can  follow  continuously  along  the 
equilibrium  path.  Nevertheless,  useful  information  can  be 
obtained  by  studying  the  static  solution  paths.1  The  lack  of  inertia 
allows  us  to  simplify  the  analysis,  since  there  is  no  need  to 


1  "Static”  refers  to  the  mode  of  recovery  of  information  on  snap-through,  not 

an  actual  physical  scenario. 


consider  the  effects  of  varying  loading  rates.  Moreover,  the 
unstable  equilibrium  path  between  the  stable  regions  discloses 
information  that  is  relevant  to  the  dynamic  case.  We  mapped  the 
static  equilibrium  path  by  using  a  pseudo-arc-length  procedure  to 
traverse  the  unstable  region. 

The  non-linear  continuum  finite  element  method  recovers 
information  that  is  unavailable  from  methods  based  on  the  beam 
or  plate  (restricted  kinematics)  equations.  In  particular,  there 
exists  a  non-trivial  deformed  configuration,  at  a  temperature 
below  the  zero-stress  reference  temperature  0ref,  that  will  not 
experience  snap-through  at  any  applied  static  load  value.  The  value 
of  this  temperature  depends  on  the  initial  beam  geometry,  the 
boundary  conditions,  and  the  pattern  of  applied  load.  Simulations 
suggest  that  any  beam  that  is  cooled  below  this  temperature  no 
longer  experiences  snap-through. 

The  subsequent  sections  in  this  paper  are  concerned  with 
developing  tests  to  establish  this  lower  bound  on  thermoelastic 
instability  and  determining  whether  this  limit,  which  is  derived 
from  purely  mechanical  considerations,  is  sufficient  to  eliminate 
dynamic  snap-through  in  beam  simulations  that  incorporate  full 
thermomechanical  coupling. 

3.  Numerical  experiments.  The  quasi-static  analysis 

The  zero-stress  reference  temperature  in  all  simulations  is 
0ref  =  300  K.  Note  that  the  adiabatic  staggered  scheme  is  imple¬ 
mented  in  terms  of  absolute  temperature,  so  0  refers  to  absolute 
temperature  and  AT  =  0—0 ref  refers  to  the  temperature  relative 
to  the  reference  temperature. 

Applied  load  versus  deflection  curves  were  extracted  for  all 
beams  at  various  temperature  variations  AT.  The  R=762  mm 
beam  (Fig.  4)  shows  relatively  little  variation  in  load-deflection 
curve  behavior  with  temperature  variation.  This  plot  is  represen¬ 
tative  for  the  results  corresponding  to  beams  with  larger  curva¬ 
tures.  Fig.  5  is  representative  for  beams  with  smaller  curvatures. 

The  key  observation  made  by  examining  Fig.  5  (and  confirmed 
for  the  other  six  beams  for  which  we  do  not  present  the  results 
here)  is  that  for  each  beam  there  exists  a  critical  temperature 
below  which  the  beam  no  longer  experiences  snap-through  at  any 
load  value  for  a  given  conservative  loading  pattern,  and  exhibits 
instead  a  monotonic  dependence  of  the  load  on  deflection.  We  will 
simply  call  this  critical  value  the  monotonic  temperature.  The  load- 
deflection  curves  for  this  particular  boundary  value  problem  also 
possess  a  center  point  where  all  the  curves  obtained  for  tempera¬ 
tures  above  the  monotonic  temperature  intersect;  this  center 
point  coincides  with  the  limit  point  where  dP/dd  =  0  on  the 
monotonic  temperature  load-deflection  curve  itself.  Results  for  Beam 
1  presented  in  Fig.  4  do  not  include  the  monotonic  curve  since  for 
this  particular  system  the  temperature  corresponding  to  it  was 
non-physical  (in  a  mathematical  sense  however,  it  does  exist). 
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Fig.  4.  Beam  1  (R=762  mm).  Pseudo-arc-length  equilibrium  paths  for  different  AT. 


The  method  presented  in  this  paper  indirectly  estimates  the 
lower  bound  under  the  assumption  of  snap-through  to  a  sym¬ 
metric  solution  branch.  This  method  relies  on  tracing  the  full 
equilibrium  path  and  on  systematic  numerical  experimentation 
to  identify  the  temperature  for  which  the  above  mentioned  path 
is  monotonic.  In  this  approach,  we  sweep  through  the  range  of 
temperatures  of  interest  and  check  the  monotonicity  of  the 
solution. 

Although  apparently  similar  results  have  been  known  for 
many  years  [5],  note  that  the  snap-free  configuration  disclosed 
in  such  studies  is  simply  the  configuration  of  an  unbuckled, 
initially  flat  beam  or  plate,  which  obviously  will  not  experience 
snap-through  under  any  magnitude  of  lateral  loading  (if  no 
temperature  variation  is  applied).  Similarly,  the  center  point  in 
these  studies  is  simply  the  unloaded  initial  configuration.  Note 
that  if  the  temperature  is  raised  above  a  critical  buckling  tem¬ 
perature,  such  a  member  will  deform  into  a  curved  shape  that  is 
itself  subject  to  snap-through  [5].  Numerous  studies  have  been 
devoted  to  mapping  this  and  the  many  subsequent  post-buckled 
states  [22,28].  Our  results  are  distinct  from  these  studies,  since 
the  underlying  problem  is  different:  here,  rather  than  an  initially 


flat  plate  or  beam  which  must  have  some  initial  thermal  defor¬ 
mation  to  experience  snap-through,  we  have  an  initially  curved 
beam  with  fixed  supports  that  displays  snap-through  at  the  zero- 
stress  reference  temperature. 

Fig.  6  indicates  that  the  monotonic  temperature  approaches 
0ref  as  fc->  0,  as  expected.  For  small  curvatures  the  lower  bound 
varies  slowly  with  the  curvature.  More  significant  variations  are 
observed  for  larger  curvatures.  For  very  large  curvatures  the 
method  presented  here  does  not  apply  since  taller  arches  will 
buckle  asymmetrically.  The  snap-through  loads  themselves  vary 
linearly  with  beam  curvature  at  AT  =  0  K,  a  fact  that  has  been 
determined  analytically  at  small  deformation  and  which  appears 
to  continue  to  hold  at  large  k\  however,  snap-through  loads  no 
longer  vary  linearly  with  temperature  away  from  the  zero-stress 
reference  temperature  (Fig.  7).  For  beams  with  very  low  initial 
curvature,  if  the  temperature  is  too  low,  the  beam  will  be  straight 
and  in  tension  and  no  snap  will  be  experienced. 

We  can  observe  several  general  trends  in  the  results. 
The  relative  influence  of  temperature  on  the  load-deflection 
behavior  of  the  beam  diminishes  as  the  curvature  k  =  1  /R  of 
the  beam  increases.  For  the  case  of  R=254mm,  the  effect  of 
even  a  temperature  change  of  AT  =  +100  1<  is  so  small  that  the 


Fig.  6.  Monotonic  temperature  varies  non-linearly  with  curvature. 


Fig.  7.  Snap-through  load  does  not  vary  linearly  with  curvature  away  from  AT  =  0  K. 
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load-deflection  curves  for  different  temperatures  would  be  diffi¬ 
cult  to  distinguish  on  a  graph;  the  snap-through  load  at  dP/dd  =  0 
barely  changed.  By  contrast,  the  influence  of  temperature  on  the 
R=  5080  mm  beam  is  significant  enough  that  AT=+100I<  is 
sufficient  to  more  than  triple  the  snap-through  load.  This  effect 
can  be  seen  in  the  decrease  in  the  magnitude  of  the  slope  with 
increasing  curvature  shown  in  plots  of  snap-through  load  versus 
temperature  in  Fig.  8.  Note  that  the  effect  of  temperature  will  be 
more  significant  for  materials  that  are  more  sensitive  to  tem¬ 
perature,  e.g.,  aluminum. 

Fig.  9  combines  these  results  and  shows  the  snap-through  load 
as  a  function  of  both  the  geometry  (radius)  and  temperature.  The 
contour  line  represents  the  snap-through  boundary,  i.e.,  the  limit 
(monotonic)  temperature  beyond  which  the  beam  does  not 
experience  snap-through.  An  alternate  method  that  can  be  used 
in  obtaining  the  snap  load  and  the  monotonic  temperature  for 
shallow  arches  is  presented  in  [29],  which  shows  a  very  good 
comparison  with  the  results  obtained  in  this  paper. 

We  also  perform  a  study  on  the  effect  of  beam  width  and 
thickness  on  the  snap-through  load  and  the  lower  bound  tem¬ 
perature.  The  study  shows  that  varying  the  beam  width  (b)  for  the 


Absolute  Temperature  [K] 


Fig.  8.  Snap-through  loads  for  all  beams  at  various  temperatures. 


same  beam  thickness  (h)  does  not  influence  the  lower  bound 
temperature;  the  snap-through  load,  however,  increases  linearly 
as  the  beam  width  increases  (Fig.  10).  As  we  vary  the  beam 
thickness,  the  snap-through  load  and  the  lower  bound  tempera¬ 
ture  are  both  higher  for  thicker  beam  (Figs.  11  and  12). 

The  pseudo-arc-length  solution  sometimes  jumps  from  the 
equilibrium  path  that  shows  the  expected  symmetric  deformation 
of  the  beam  (Fig.  13)  to  another  path  where  the  beam  deforms 
asymmetrically  (Fig.  14).  This  asymmetric  buckling  mode  has 
been  previously  discussed  [5].  While  the  second  buckling  mode  is 
negligible  if  the  initial  curvature  is  small  enough  [27],  this  mode 
clearly  cannot  be  neglected  for  large  initial  beam  curvatures. 

Nevertheless,  a  unique  lower  bound  on  instability  for  a  given 
beam  and  load  pattern  can  still  be  established  if  the  second  mode 
only  occurs  above  the  monotonic  temperature.  If  the  entire 
unstable  region  itself  is  eliminated  from  the  static  solution  path 
at  the  monotonic  temperature  -  including  the  higher  buckling 
modes  -  then  the  lower  bound  on  instability  is  unique  for  a  given 
load  pattern  and  initial  member  shape.  No  static  simulation  that 
we  conducted  contradicted  this  conjecture,  and  dynamic  simula¬ 
tion  results  indicated  that  higher  modes  tend  to  be  activated  at 
higher  temperatures  but  disappear  as  temperature  is  lowered, 
suggesting  that  the  lower  bound  is  in  fact  unique. 
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3.1.  The  monotonic  temperature  and  the  energy  derivative  tests 

The  coincidence  of  the  limit  point  of  the  load-deflection  curve 
at  the  monotonic  temperature  with  the  center  point  of  the  graph 
suggests  that  there  is  something  special  about  this  point,  and 
inspection  of  the  graphs  suggests  that  this  point  is  likely  the 
inflection  point  of  the  curve  (Figs.  15  and  16).  Numerical  evalua¬ 
tion  of  the  derivative  c^P/Sd2  confirms  this.  Note  that  these 
results  are  obtained  through  numerical  approximation,  therefore 
the  coincidence  points  are  not  exact.  However,  we  prove  through 
analytical  studies  that  this  hypothesis  holds  for  a  simple  case  [26]. 
In  the  case  of  the  symmetric  point-loaded  beam,  we  can  construct 
a  straightforward  physical  explanation  of  the  existence  of  this 
center  point  that  will  help  guide  our  understanding  of  more 
general  load  cases  and  geometries. 

If  the  beam  stays  symmetric  about  the  midpoint,  the  slope  of 
the  beam  at  the  midpoint  must  always  be  horizontal,  and  by 
equilibrium  the  axial  internal  force  at  the  midpoint  must  be  equal 
to  the  horizontal  reaction  force  at  the  supports.  We  intuitively 
expect  that  the  unstable  configuration  with  the  maximum  hor¬ 
izontal  reaction  force  is  the  maximally  unstable  configuration, 
that  is,  the  unstable  configuration  that  would  move  to  a  stable 
configuration  with  the  maximum  possible  kinetic  energy  relative  to 
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all  other  unstable  configurations  (if  the  beam  were  initially  held 
perfectly  fixed  at  this  configuration  and  then  released). 

The  maximum  horizontal  reaction  does  coincide  with  the 
center  point  of  the  graphs  (Fig.  17),  indicating  that  the  center 
point  corresponds  to  the  maximally  unstable  configuration.  This 
does  not  mean  that  the  horizontal  reaction  must  be  less  than  or 
equal  to  zero  at  the  monotonic  temperature  in  order  to  avoid 
instability,  as  one  might  assume:  Fig.  17b  shows  that  the  max¬ 
imum  horizontal  reaction  can  still  be  compressive  at  the  mono¬ 
tonic  temperature. 

The  center  point  is  a  function  of  the  specific  boundary  value 
problem  that  we  chose  to  explore.  We  would  like  to  develop  a 
more  general  test  to  determine  the  maximally  unstable  config¬ 
uration  for  problems  where  the  load-deflection  curve  may  not  be 
so  well-behaved.  The  following  heuristic  argument  will  point  us 
toward  such  a  test.  For  this  particular  problem,  the  external  work 
can  be  computed  as  the  integral  of  the  product  of  the  concen¬ 
trated  load  P(d)  and  the  deflection  at  the  midpoint  d. 

The  second  derivative  of  energy  gives  the  standard  static 
stability  test,  while  the  third  derivative  of  energy  gives  new 
information  that  should  locate  the  center  point.  The  third  deri¬ 
vative  of  energy  allows  us  to  suggest  a  mathematical  definition  of 
the  maximally  unstable  configuration  and  the  monotonic  tem¬ 
perature.  Since  the  value  of  cPE/dd2  can  be  considered  a  measure 
of  the  degree  of  instability  present  in  the  system,  we  would  expect 
that  a  point  where  its  derivative  cPE/dd3  is  zero  would  be  an 
instability  extremum  that  is  local  to  the  unstable  region.  The 
maximally  unstable  configuration  would  be  the  instability  extre¬ 
mum  in  the  range  where  cPE/dd2  <  0.  If  we  are  able  to  shrink  this 
range  by  lowering  the  temperature,  then  in  the  limit  where  the 
unstable  domain  shrinks  to  zero  we  have  both  cPE/dd3  =  0  and 
82E/dd2  =  0  at  the  same  configuration.  The  temperature  that 
generates  this  configuration  and  satisfies  these  conditions  is  the 
monotonic  temperature.  In  other  words,  the  temperature  that 
shrinks  the  unstable  region  to  a  point,  for  which  the  maximally 
unstable  ( 83E/8d 3  =0)  and  the  minimally  stable  ( d2E/dd 2  =0)  con¬ 
figurations  are  necessarily  one  and  the  same,  is  the  monotonic 
temperature.  Note  that  this  is  similar  to  the  use  of  higher 
derivatives  in  continuation  methods  [22],  though  some  simplifi¬ 
cations  of  these  more  general  (and  computationally  expensive) 
techniques  may  be  possible  for  thermomechanical  problems,  as 
discussed  below. 

The  energy  Emt  is  obtained  from  the  finite  element  code,  and 
its  numerical  derivatives  can  be  approximated  by  simple  differ¬ 
ence  formulas.  We  use  this  data  to  test  the  above  hypothesis  and 


we  do  find  that  the  third  numerical  derivatives  of  energy  are  zero  at 
the  center  point,  and  the  second  derivative  of  energy  is  zero  at  the 
same  point  for  the  monotonic  temperature  (Fig.  15).  This  observa¬ 
tion  holds  for  beams  of  significantly  different  curvatures  (Fig.  16). 

For  the  case  of  the  distributed  load  (Fig.  2b),  a  global  center 
point  no  longer  exists  (Fig.  18).  The  curves  for  lower  temperatures 
possess  a  center  point  but  the  curves  at  elevated  temperatures  do 
not  pass  through  it.  Similar  to  the  case  of  the  concentrated  load, 
the  curved  beam  acts  like  an  arch  under  distributed  load,  resisting 
primarily  through  axial  stiffness  rather  than  through  weaker 
flexural  stiffness.  However,  this  “arch”  is  not  geometrically 
perfect  and  suffers  from  local  flexural  buckling  (Fig.  20a-c)  at 
load  values  below  those  achieved  in  the  case  of  the  concentrated 
load. 

These  higher  buckling  modes  may  deflect  upward  at  the 
midpoint  even  as  global  stability  is  lost,  so  we  can  no  longer 
assume  that  dP/dd  =  0  indicates  a  snap-through  point,  at  least  at 
elevated  temperatures.  Similarly,  d3E/8d3  =  0  does  not  occur  in 
the  distributed  load  case  at  high  temperatures  and  can  therefore 
no  longer  be  related  to  the  maximally  unstable  configuration  as  it 
can  in  the  concentrated  load  case.  Nonetheless,  the  energy 
derivatives  are  still  zero  at  the  monotonic  temperature  (Fig.  19). 


Total  Force  Applied  [N] 

Fig.  18.  Beam  5  (R= 2438.4  mm).  Load-deflection  curves  for  distributed  loads:  no 
center  point  at  elevated  temperatures. 


a  b 


Fig.  17.  Beam  8  (R— 5080  mm).  Maximum  horizontal  reaction  occurs  at  center  point:  (a)  load-deflection  diagram  and  (b)  beam  horizontal  reaction. 
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The  distributed  load  case  reinforces  the  earlier  observation 
that  the  value  of  the  monotonic  temperature  depends  on  the  pattern 
of  applied  load:  for  the  R=2438.4  mm  beam,  the  symmetric  point 
load  leads  to  a  monotonic  temperature  of  A T  =  -22  K,  while  for 
the  distributed  load  case  the  monotonic  temperature  is  approxi¬ 
mately  AT  =  -28  K.  The  decreased  monotonic  temperature  for  the 
distributed  load  case  indicates  that  it  is  possible  to  find  a  loading 
such  that  a  “worst-case"  monotonic  temperature  is  obtained  and 
a  lower  bound  on  the  thermoelastic  instability  for  all  possible 
loadings  of  a  system  is  established.  This  requires  the  solution 
of  an  inverse  problem  to  determine  the  worst-case  loading, 
which  entails  more  complexity  than  is  typically  worth  the  effort, 


Fig.  19.  Beam  5  (R=2438.4  mm).  Energy  derivatives  are  still  zero  at  monotonic 
temperature  despite  the  absence  of  a  global  center  point. 


especially  if  the  actual  load  patterns  a  structure  is  subject  to  can 
be  known  in  advance  with  reasonable  certainty.  This  concept  is 
not  developed  further  in  this  paper,  and  all  monotonic  tempera¬ 
tures  reported  are  limited  to  the  specific  load  cases  considered. 

In  all  these  examples,  due  to  symmetry,  the  derivatives  with 
respect  to  d  (the  midpoint  deflection)  provide  clear  information 
about  the  stability  of  the  system.  In  general,  there  will  not  be  a 
simple  line  of  symmetry  that  will  allow  us  to  obtain  meaningful 
information,  as  is  apparent  from  the  higher  buckling  mode  load- 
deflection  curves  in  Fig.  20.  It  would  be  preferable  to  develop  an 
analytical  test  for  the  third  derivative  of  energy  that  will  let  us 
circumvent  the  process  of  taking  numerical  derivatives  entirely 
and  allow  us  to  determine  the  monotonic  temperature  for  any 
given  system. 

Before  assuming  that  the  monotonic  temperature  obtained 
from  static  simulations  is  truly  effective  in  eliminating  snap- 
induced  vibrations,  we  must  first  demonstrate  that  it  eliminates 
snap-through  in  dynamic  systems.  The  next  section  is  devoted  to 
establishing  the  validity  of  generalizing  from  the  static,  mechan¬ 
ical-only  solution  to  the  dynamic,  thermomechanically  coupled 
case.  Then,  we  outline  a  test  for  establishing  the  monotonic 
temperature  for  arbitrary  loading  and  member  geometries  in 
the  final  section. 


4.  Dynamic  simulation  of  snap-through 

In  order  to  determine  if  the  monotonic  temperature  estimated 
via  static  simulations  successfully  eliminates  snap-through, 
dynamic  simulations  are  conducted.  We  present  here  the  results 
obtained  for  one  beam  only  (R=  1828.8  mm),  with  the  same 
reference  temperature  0ref  =  300  K  used  in  the  static  simulations. 
Consequently,  the  beam  number  is  no  longer  mentioned  in  the 
figure  captions.  This  particular  beam  has  a  large  initial  curvature 
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Fig.  20.  Beam  5  (R— 2433.4  mm)  with  distributed  load:  (a)  AT  =  0  K,  (b)  AT  =  5  I<  and  (c)  AT  =  10  K. 
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but  still  possesses  an  estimated  monotonic  temperature  that 
could  be  easily  obtained  in  the  laboratory  (AT  =  -38  K  relative 
or  0  =  -11  °C/12  °F  absolute).  A  concentrated  force  is  applied  at 
an  effectively  quasi-static  rate  at  midspan  and  then  held  constant 
at  a  fixed  target  value  greater  than  the  static  snap-through  load  Ps. 
The  target  load  is  P=2.0N  for  temperatures  below  T=30K,  but 
the  increase  in  snap-through  load  at  increased  temperature 
necessitated  higher  target  loads  of  P=  2.5  and  3.0  N  for  T=  50 
and  100  K  respectively  (Fig.  21). 

Initial  static  simulations  were  conducted  to  obtain  the  ther¬ 
mally  deformed  initial  configuration  of  the  beam  at  each  target 
temperature  prior  to  dynamic  simulation  (Fig.  22).  The  initial 
temperature  was  applied  as  a  perfectly  uniform  nodal  initial 
condition  and  the  surfaces  of  the  beam  were  treated  as  perfectly 
insulated.  Though  this  zero-heat-flux  boundary  condition  is  not 
feasible  in  an  experimental  setting,  it  allows  us  to  isolate  the 
temperature  change  due  to  thermomechanical  coupling  in  our 
simulations.  Purely  mechanical  dynamic  simulations  [20]  had 
earlier  indicated  that  a  problem-specific  critical  time  step  of 
At=10~4s  or  smaller  must  be  used  for  accurate  modeling  of 


Fig.  21.  Ramp  load  target  values. 


snap-through.  Values  of  At  above  this  critical  value  could  lead  to 
spurious  solutions  that  are  clearly  non-physical  but  nonetheless 
numerically  converged. 

The  monotonic  temperature  estimated  from  static  simulations 
successfully  eliminated  snap-through  from  dynamic  simulations 
at  quasi-static  load  rates.  The  maximum  amplitude  of  displace¬ 
ment  does  clearly  approach  zero  as  the  temperature  approaches 
the  monotonic  temperature,  and  effectively  reaches  zero  for 
AT  =  -50  K  (Fig.  23).  The  maximum  displacement  approaches  a 
limit  as  the  temperature  increases;  however,  the  maximum 
kinetic  energy  over  all  time  values  continues  to  increase  linearly 
as  the  temperature  increases  (Fig.  24).  This  “excess”  kinetic 
energy  is  present  in  asymmetric  deformation  modes  not  apparent 
from  examination  of  the  midpoint  displacement  alone. 

Note  that  the  ramp  load  continues  to  increase  prior  to  reaching 
a  constant  value  after  snap-through,  hence,  the  centerline  of  the 
vibrations  moves  downward  on  the  plot  prior  to  f=2.0  s  (Fig.  26). 


Fig.  23.  Maximum  amplitude  of  vibrations  at  various  temperatures.  Vibrations 
disappear  near  monotonic  temperature  and  level  off  at  elevated  temperature. 


Fig.  24.  Maximum  kinetic  energy  at  various  temperatures.  Kinetic  energy  disap¬ 
pears  near  monotonic  temperature  and  increases  linearly  at  elevated  temperatures. 
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This  component  of  the  vibration  must  be  filtered  out  in  order 
to  determine  the  amplitudes  of  vibration  correctly.  This  was 
done  by  fitting  a  sixth-order  polynomial  to  the  post-snap  mid¬ 
point  displacement  time  history  data  and  finding  the  value 
of  the  maximum  difference  between  the  smooth  curve  and  the 
displacement. 

The  y-direction  displacement  at  the  midpoint  is  represented  in 
Figs.  25-27  as  a  function  of  time  for  various  initial  temperatures 
and  the  static  simulation  results  are  superposed  at  the  same  scale. 
Thermoelastic  damping  is  clearly  apparent  in  these  plots;  by 
comparison,  purely  mechanical  simulations  of  snap-through 
using  the  trapezoidal  Newmark  method  showed  no  damping  of 
the  post-snap  vibrations,  as  is  expected  for  this  energy-conserving 
time  integration  algorithm.  Visual  inspection  of  the  graphs  also 
indicates  that  the  large-amplitude  vibrations  characteristic  of 
snap-through  are  strongly  attenuated  below  the  monotonic 
temperature  of  AT  =  -38  K  (Fig.  26),  dropping  to  nearly  zero  at 
AT  =  —50  K  (Fig.  27). 

A  measure  of  the  asymmetry  of  vibration  can  be  constructed 
by  plotting  the  difference  between  the  y-direction  displacement 


Time  [s] 


Fig.  25.  Displacement  at  AT  =  0  K. 


Time  [s] 

Fig.  27.  Displacement  at  AT  =  —50  K. 


w(x)  at  the  points  x  =  L/ 4  and  x  =  3T/4.  This  difference  will  be 
zero  if  the  vibration  is  symmetric,  or  otherwise  can  give  an 
estimate  of  the  frequency  and  magnitude  of  the  asymmetric 
vibration.  For  AT  =  0  K  and  below,  the  post-snap  vibration 
remains  perfectly  symmetric  for  the  beam  considered.  For 
AT  =  30  K,  the  post-snap  vibration  is  initially  symmetric,  but  after 
the  load  reaches  a  constant  value  the  response  develops  an 
antisymmetric  vibration  (Fig.  28).  This  vibration  damps  out 
relatively  quickly  at  AT  =  30  K,  but  at  AT  =  50  K  a  similar  mode 
of  vibration  is  present  at  a  higher  amplitude  and  for  a  longer 
duration  (Fig.  29).  For  AT  =  100  K,  the  beam  buckles  asymmetri¬ 
cally  immediately  prior  to  the  static  snap-through  load  and 
remains  in  a  fully  asymmetric  mode  after  snap-through  (Fig.  30). 

Qualitatively,  the  AT  =  30  K  and  AT  =  50  K  cases  combine  a 
symmetric  mode  (Fig.  31a)  and  an  extremely  low  frequency 
asymmetric  mode  (Fig.  31b).  The  AT  =  100  1<  case  buckles  asym¬ 
metrically  prior  to  snap-through  (Fig.  32a)  and  remains  in  a  high- 
frequency  asymmetric  mode  that  oscillates  at  an  angle  to  the 
vertical  following  buckling  (Fig.  32b).  As  this  asymmetric  mode  is 
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Fig.  29.  Low-frequency  asymmetric  vibration  sustains  longer  for  AT  =  50  K. 


a 


Fig.  31.  Mode  shapes  for  AT  =  30  and  50  K:  (a)  symmetric  and  (b)  asymmetric. 

damped,  a  new  mode  emerges  (following  t=3.4  s)  similar  to  the 
low-frequency  asymmetric  mode  encountered  at  AT  =  30  K  and 
AT  =  50  K. 


Fig.  32.  AT=100I<:  (a)  initial  asymmetric  buckling  and  (b)  asymmetric  "side¬ 
ways”  post-buckling  mode. 


b 


Fig.  33.  Modes  at  unloading  for  AT  =  0  K:  (a)  symmetric  and  (b)  traveling  wave. 


Unloading  the  beam  from  an  initially  stationary  deformed 
configuration  reveals  that  the  vibration  about  the  initial  unde¬ 
formed  configuration  is  of  significantly  lower  amplitude  than  the 
vibration  about  the  deformed  configuration.  Although  the  initial 
post-snap  unloading  response  begins  as  a  symmetric  vibration,  it 
develops  a  small  asymmetry  that  quickly  grows  in  magnitude. 
While  the  initial  transient  mode  shape  is  essentially  similar  to  the 
low-order  symmetric  mode  of  the  loading  response  (Fig.  33a), 
inspection  of  the  deformed  configurations  reveals  that  the  asym¬ 
metric  mode  behaves  like  a  low-frequency  "wave”  reflecting  back 
and  forth  between  the  supports  (Fig.  33b).  Fig.  34  shows  the 
loading  and  unloading  responses. 

These  various  modes  interact  non-linearly  in  the  finite-defor¬ 
mation  model,  and  the  resulting  transient  behavior  can  become 
difficult  to  quantify,  especially  as  higher  modes  are  activated  and 
the  post-snap  motion  increases  in  complexity.  However,  the 
higher  modes  diminish  with  decreasing  temperature  and  even¬ 
tually  disappear  at  the  monotonic  temperature,  making  the 
monotonic  temperature  a  solid  point  of  reference  in  an  otherwise 
convoluted  post-buckling  regime. 

The  monotonic  temperature  also  holds  for  the  distributed  load 
case.  The  R=2438.4mm  beam  displays  the  same  sort  of  asym¬ 
metric  buckling  prior  to  the  static  snap-through  load  at  AT  =  0  K 
under  distributed  loading  that  is  seen  at  AT=100K  for  the 
R=  1828.8  mm  beam  under  concentrated  loading.  At  the  esti¬ 
mated  monotonic  temperature  nonetheless,  post-snap  vibrations 
are  attenuated  as  effectively  in  the  distributed  load  case  as  they 
are  in  the  concentrated  load  case  (Fig.  35). 

There  are  some  caveats  on  the  generality  of  the  monotonic 
temperature.  The  softening  response  seen  in  the  static  solution 
paths  even  for  temperatures  below  the  monotonic  temperature 
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Fig.  34.  Beam  3  (R=  1828.8  mm).  Asymmetry  of  loading  and  unloading  response 
at  AT  =  0  K. 


Time  [s] 

Fig.  35.  Beam  5  (R=2438.4  mm).  Static  monotonic  temperature  estimate  holds  for 
distributed  load  case. 

can  still  lead  to  snap-like  behavior  if  the  loading  rate  is  high 
enough  to  force  the  system  to  dynamically  jump  off  the  static 
solution  curve  as  the  system  softens  with  increasing  load.  More 
generally,  the  system  must  be  dissipative  for  the  temperature  to 
remain  below  the  monotonic  temperature.  This  is  true  in  our  case, 
but  it  is  no  longer  valid  if  body  heat  sources  or  applied  surface 
heat  fluxes  are  present.  We  are  therefore  limited  to  saying  that 
the  monotonic  temperature  provides  a  lower  bound  below  which 
snap-through  will  never  occur  at  any  applied  load  level  if  the 
loading  rate  is  effectively  quasi-static  and  the  system  is 
dissipative. 

The  monotonic  temperature  estimated  from  static  simulations 
is  not  perfectly  precise;  it  does  not  completely  eliminate  the 
oscillations  even  for  dynamic  systems  with  effectively  quasi¬ 
static  loading.  Some  of  these  errors  are  due  to  the  simple 
approach  that  was  undertaken  to  find  this  temperature,  where 
the  monotonicity  of  the  curve  was  determined  only  approxi¬ 
mately.  Another  source  for  the  discrepancy  is  due  to  the  physical 
fact  that  the  ramp  loading  can  never  be  truly  quasi-static  in  a 
dynamic  simulation.  The  amplitude  of  vibration  will  go  to  zero 
nearer  the  monotonic  temperature  (refer  to  Fig.  23)  if  the  loading 
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rate  is  decreased,  but  there  will  always  exist  some  small  dynamic 
jump  as  the  system  approaches  the  limit  point  on  the  static  load- 
deflection  curves. 

One  additional  practical  observation  can  be  made.  We  have 
modeled  the  curved  beam  with  an  edge  that  is  free  to  move  in  the 
z-direction  (out-of-plane)  between  the  fixed  supports.  In  many 
practical  engineering  applications,  such  as  stressed-skin  mono- 
coque  aircraft  construction,  the  z-direction  will  be  restrained  as 
well,  and  we  expect  this  additional  boundary  condition  to  have  a 
stabilizing  effect.  These  additional  boundary  conditions  increase 
the  snap-through  load  (Fig.  36)  as  well  as  the  monotonic  tem¬ 
perature  itself  (Fig.  37).  The  monotonic  temperatures  obtained 
from  the  static  simulations  without  these  boundary  conditions 
are  therefore  conservative  values  relative  to  systems  that  include 
such  constraints. 

Note  that,  in  the  static  simulations,  we  have  omitted  heat 
conduction  and  treated  temperature  as  a  purely  mechanical 
expansion  or  contraction  of  the  material,  which  provided  a 
valid  estimate  of  the  snap-free  monotonic  temperature  for  the 
dynamic  thermomechanically  coupled  case.  This  analogy  between 
mechanical  strain  and  thermal  strain  indicates  that  an  applied 
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Fig.  36.  Beam  3  (R= 1828.8  mm).  Full  z-direction  edge  boundary  condition 
increases  snap-through  load. 


Fig.  37.  Beam  3  (R= 1828.8  mm).  Full  z-direction  edge  boundary  condition 
increases  monotonic  temperature. 
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mechanical  prestress  will  have  a  mechanical  effect  similar  to 
lowering  the  temperature.  For  example,  for  the  unloaded 
R=  1270  mm  beam,  the  reaction  at  the  supports  has  no  vertical 
component  and  has  a  horizontal  component  of  Px  =  9.42  N  at  the 
monotonic  temperature  AT  =  —82  K;  the  horizontal  reaction  is 
PX=15.4N  for  Ar=-100K.  Given  the  out-of-plane  depth  of 
d=28.7  mm  and  beam  thickness  t=1.27mm,  and  assuming  a 
uniform  stress  over  the  cross-section,  the  axial  stress  in  the  beam 
is  therefore  0.4225  MPa  for  AT  =  -100  K,  which  is  only  a  minute 
fraction  of  the  yield  stress  of  <jy  =  250  MPa  for  A36  structural 
steel.  Although  AT  =  -100  K  is  clearly  an  unrealistic  temperature, 
a  =  0.4225  MPa  is  a  perfectly  realistic  prestress.  Consequently,  an 
applied  mechanical  prestress  will  have  a  mechanical  effect  similar 
to  lowering  the  temperature,  and  therefore  we  could  eliminate 
quasi-static  snap-through  by  applying  a  purely  mechanical  pre¬ 
tension  sufficient  to  achieve  a  tensile  strain  identical  to  or  greater 
than  that  produced  by  the  monotonic  temperature.  Moreover,  it  is 
also  possible  to  apply  additional  initial  strain  sufficient  to  cancel 
the  effects  of  thermal  expansion,  so  that  a  heated  beam  can  be 
“snap-proofed”  for  temperatures  below  a  given  elevated  target 
temperature. 

In  practical  aerospace  applications,  significant  prestress  is 
introduced  to  thin  structural  members  to  provide  increased 
structural  stiffness,  independent  of  any  thermal  considerations. 
This  prestress  is  sufficient  to  introduce  tensile  strains  equivalent 
or  greater  than  those  introduced  by  the  monotonic  temperature, 
as  is  apparent  from  the  above  example.  It  is  therefore  unlikely 
that  typical  aerospace  structures  will  experience  snap-through 
under  standard  operating  conditions.  However,  as  the  operating 
conditions  of  the  structure  become  more  severe,  as  in  hypersonic 
flight,  much  larger  thermal  loads  can  be  expected.  Additional 
mechanical  pre-tension  could  therefore  prove  useful  in  providing 
additional  protection  against  snap-through  in  aerospace  struc¬ 
tures  in  extreme  operating  environments,  although  more  research 
would  be  required  to  determine  the  effectiveness,  applicability, 
and  limitations  of  this  approach. 


5.  Maximum  instability  and  bounds  on  snap-through 

In  the  previous  sections,  we  have  established  that  (1)  the 
monotonic  temperature  represents  a  lower  bound  on  static 
instability  in  curved  beams  subject  to  snap-through,  and  (2)  for 
the  simple  symmetric  beam  with  a  concentrated  load,  the  second 
and  third  derivatives  of  energy  with  respect  to  the  midpoint 
displacement  are  both  zero  at  the  snap-through  point  at  this 
temperature.  We  surmised  that  the  unstable  region  shrinks  as  the 
temperature  is  lowered  until  it  becomes  a  point  at  the  monotonic 
temperature,  and  that  at  this  point  the  minimally  stable  config¬ 
uration  (52£/5d2  =  0)  and  maximally  unstable  configuration 
(i tPE/dd 3  =  0)  are  identical.  This  shrinkage  of  the  unstable  region 
to  a  point  as  temperature  is  lowered  is  illustrated  conceptually  in 
Fig.  38.  We  also  observed  that,  although  multiple  modes  may 
exist  corresponding  in  some  cases  to  multiple  paths  within  the 
unstable  region  (Fig.  39),  in  static  simulations  these  modes  tend 
to  disappear  as  the  monotonic  temperature  is  approached. 

These  observations  suggest  a  mathematical  definition  of  the 
monotonic  temperature.  We  do  not  attempt  to  prove  the  follow¬ 
ing  conjectures,  but  the  simulation  results  support  them.  If  we  are 
able  to  shrink  the  unstable  region  along  the  equilibrium  path  of 
our  system  by  altering  some  parameters  (such  as  temperature), 
then  the  configuration  where  the  critical  parameters  are  such  that 
snap-through  no  longer  occurs  is  exactly  the  point  where  the 
minimally  stable  configurations  at  the  boundaries  of  the  unstable 
region  coincide  with  the  maximally  unstable  configuration.  Such 
conditions  should  therefore  give  us  a  set  of  equations  that  provide 
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a  bound  of  the  domain  of  unstable  behavior.  Without  attempting 
to  prove  them,  we  state  here  these  conditions. 

Conjecture.  Let  d  be  the  displacement  field  in  a  continuum,  and  let 
p  be  a  vector  of  bifurcation  parameters,  which  control  the  size  of  the 
unstable  region.  Assume  that  the  system  under  consideration  lias 
only  one  unstable  region.  Let  C,  A,  B  be  the  first,  second  and  third 
derivatives  of  the  energy.  The  unstable  region  ceases  to  exist  for  the 
particular  configuration  (d*,p*)  at  which,  for  all  virtual  displacement 
fields  Sd,  the  functionals  G,  A  and  B  are  null :  G(d*,p*,Sd)  =  0, 
A(d*,p*,Sd)  =  0,  and  B(d*,p*,Sd)  =  0. 

If  we  use  a  finite  element  approximation  of  the  functionals 
from  a  continuum  problem,  we  obtain  the  version  of  the  above 
conditions  in  matrix  form:  G(d*,p*)  =  0,  detA(d*,p*)  =  0,  and 
det  B(d*,p*)  =  0.  In  this  discretized  form,  G  is  a  vector  (the 
residual  vector  in  finite  element  methods),  A  is  a  second-order 
tensor  (the  stiffness  matrix),  B  is  a  third-order  tensor,  and  d*  is  a 
vector  (of  nodal  displacements). 

A  simple  non-linear  Timoshenko  beam  model  demonstrates  the 
utility  of  these  conjectures  and  also  demonstrates  why  limited- 
deflection  models,  such  as  the  von  Karmann  equations,  are  a  priori 
incapable  of  determining  the  monotonic  temperature.  The  inter¬ 
ested  reader  is  referred  to  [26]  for  the  detailed  derivation  corre¬ 
sponding  to  the  beam  model.  Unlike  the  fully  non-linear  equilibrium 
equations,  the  third  variational  derivative  of  the  linearized  beam 
equations  is  identically  zero.  In  taking  the  Taylor  series  approxima¬ 
tion  of  the  fully  non-linear  Timoshenko  beam  equilibrium  equations 
we  have  lost  information,  and,  as  discussed  in  [26],  the  information 
we  have  lost  from  the  third  derivative  of  energy  is  exactly  the 
information  we  require  to  be  able  to  pinpoint  the  non-trivial  snap- 
free  configuration  that  exists  at  the  monotonic  temperature.  The  von 
Karmann  equations  therefore  do  not  disclose  the  existence  of  the 
non-trivial  solution  precisely  because  the  information  that  is  needed 
to  determine  it  has  been  omitted  a  priori. 
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6.  Conclusions 

We  demonstrated  the  existence  of  a  non-trivial  curved  beam 
configuration  and  a  corresponding  critical  temperature,  lowered 
with  respect  to  the  reference  temperature,  that  eliminates  the 
occurrence  of  dynamic  snap-through  under  quasi-static  loading 
rates.  This  monotonic  temperature  can  be  estimated  from  static 
solution  paths  with  reasonable  accuracy,  even  for  beams  of  large 
initial  curvature.  An  analytical  test  based  on  the  third  variational 
derivative  of  the  energy  that  could  be  useful  in  determining  the 
monotonic  temperature  for  members  of  general  geometry  and 
load  patterns  is  also  suggested.  The  failure  of  limited-deformation 
thermal  beam  and  plate  models  such  as  the  von  Karmann  equations 
to  disclose  the  existence  of  a  non-trivial  snap-free  solution  is  shown 
to  be  a  consequence  of  the  omission  of  information  about  the  third 
variational  derivative  of  energy  from  these  lower-order  theories. 

The  analogy  between  purely  mechanical  contraction  and 
thermal  strain  that  holds  for  these  numerical  experiments  sug¬ 
gests  that  an  applied  mechanical  prestress  that  introduces  an 
initial  strain  equivalent  to  the  strain  induced  by  the  monotonic 
temperature  eliminates  the  possibility  of  quasi-static  snap- 
through  in  beams  of  large  initial  curvature.  Increasing  the  initial 
prestress  in  thin-walled  structures  using  values  determined  from 
this  procedure  could  therefore  provide  a  simple  method  of 
protecting  against  large  amplitude  snap-through  vibrations 
encountered  in  by  aircraft  and  spacecraft  in  extreme  operating 
environments.  Although  this  method  would  introduce  additional 
stresses  into  the  supporting  members  and  therefore  has  limita¬ 
tions  depending  on  the  requirements  of  the  structure,  it  is  simpler 
than  other  methods  that  have  been  proposed  in  the  literature 
[30].  Further  research  would  be  needed  to  determine  the  poten¬ 
tial  effectiveness  and  practical  limitations  of  this  proposed 
method  of  alleviating  dynamic  snap-through. 
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Aircraft  structures  operating  in  severe  environments  may  experience  snap-through,  causing  the 
curvature  on  part  or  all  of  the  structure  to  invert  inducing  fatigue  damage.  This  paper  examines  the 
performance  of  beam  and  continuum  nonlinear  finite  element  formulations  in  conjunction  with  several 
popular  implicit  time  stepping  algorithms  to  assess  the  accuracy  and  stability  of  numerical  simulations  of 
snap-through  events.  Limitations  of  the  structural  elements  are  identified  and  we  provide  examples  of 
interaction  between  spatial  and  temporal  discretizations  that  affect  the  robustness  of  the  overall  scheme 
and  impose  strict  limits  on  the  size  of  the  time  step.  These  limitations  need  to  be  addressed  in  future 
works  in  order  to  develop  accurate,  robust  and  efficient  simulation  methods  for  response  prediction  of 
structures  encountering  extreme  environments. 

©  201 1  Elsevier  Ltd.  All  rights  reserved. 


1.  Introduction 

Curved  beams  or  panels  can  often  be  found  as  components  of 
complex  structures  in  civil,  mechanical,  and  aerospace  applications. 
They  may  experience  snap-through  causing  the  curvature  on  part  or 
all  of  the  structure  to  invert  due  to  a  large  inward  loading  (Fig.  1(a)) 
caused  by  extreme  loading  conditions  [1],  When  snap-through 
happens,  the  load-deflection  diagram  presents  a  jump,  e.g.,  from 
point  B  to  point  C,  as  shown  in  Fig.  1  ( b).  In  some  cases,  when  the  load 
is  reduced,  hysteresis  is  observed:  the  structure  will  snap-back  to  its 
original  configuration  but  on  a  different  path  (represented  by  the 
jump  from  point  D  to  point  A  in  the  same  figure). 

Snap-through  is  characterized  by  large  nonlinear  deformations, 
changes  in  the  system  stability  and  large  stress  reversals,  which 
accelerate  fatigue  damage.  Since  no  analytical  solutions  for  general 
systems  with  snap-through  exist,  numerical  models  that  can  capture 
these  phenomena  are  needed  in  order  to  predict  the  fatigue  life  of  the 
structure.  Among  the  numerical  techniques,  the  Finite  Element 
Method  (FEM)  provides  the  most  generality  and  can  be  applied  to 
systems  with  arbitrarily  complex  geometries.  This  paper  analyzes  the 
performance  of  several  finite  element  formulations  (2D  and  3D) 
and  the  stability  of  the  time-stepping  schemes  in  simulating  a 
curved  beam  undergoing  snap-through:  we  identify  the  important 
features  that  affect  the  numerical  accuracy  and  robustness  and  the 
region  where  the  schemes  are  stable  for  such  simulations. 
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The  snap-through  of  curved  beams  or  shallow  arches  has  been 
studied  analytically  for  particular  cases  by  numerous  authors, 
including  Murphy  et  al.  [2],  Virgin  [3],  Bradford  et  al.  [4],  and  Plaut 
and  Virgin  [5].  Snap-through  can  also  be  studied  qualitatively 
using  a  truss  system  [6].  Simplifying  assumptions  are  usually 
needed  in  finding  the  analytical  estimates.  While  we  are  paying 
the  price  of  the  lack  of  physical  intuition  by  using  the  complex 
FEM,  we  often  need  to  do  so  in  order  to  use  more  general  analysis 
techniques  that  allow  for  complex  effects  to  be  captured.  There¬ 
fore,  we  test  here  the  performance  of  this  method  when  applied 
to  a  simple  curved  beam  structure  (geometry  and  properties 
introduced  in  Section  2).  This  structure  is  representative  since  it 
exhibits  the  type  of  phenomena  we  wish  to  study  (snap-through) 
but  simple  enough  to  allow  for  comparisons  with  solutions 
obtained  with  analytical  methods  as  well. 

In  this  paper,  the  numerical  simulations  are  performed  with  the 
Finite  Element  Analysis  Program  (FEAP),  a  research  code  that 
includes  most  commonly  used  finite  elements  and  solvers  and 
provides  a  reliable  framework  for  developing  and  implementing 
new  user  formulations  [7],  The  factors  taken  into  consideration 
in  this  study  are  (1)  the  time-stepping  algorithm,  (2)  the  element 
type,  and  (3)  the  size  of  the  time  step.  The  element  formulations 
and  the  time-stepping  algorithms  discussed  are  either  available 
from  the  standard  FEAP  distribution  or  implemented  as  user 
subroutines.  The  main  criterion  used  to  verify  the  robustness  of 
the  results  is  the  energy  conservation  throughout  the  simulation. 

For  discretization,  2D  straight  beam  elements  and  3D  solid 
linear  and  quadratic  elements  are  used.  All  formulations  account 
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Fig.  1.  Snap-through  buckling  of  shallow  structures:  (a)  initial  and  post-snap 
configuration,  (b)  load-deflection  curve. 


for  large  deformations.  The  2D  beam  elements  utilized  are  (1) 
beam  elements  without  shear  deformation  with  large  displacement 
and  small  rotation  (2nd  order  theory)  with  cubic  interpolation,  and 
(2)  beam  elements  with  shear  deformation  formulation  and  large 
displacement  and  large  rotation  that  can  consider  inelastic  behavior 
for  bending  and  axial  effects  but  retain  linear  elastic  response  in 
the  transverse  shear  terms  [7].  The  3D  solid  elements  used  are 
linear  (8-nodes)  elements  with  displacement,  mixed  (B-bar),  and 
enhanced  formulations  and  quadratic  (27-nodes)  elements  with 
displacement  formulation.  Unlike  the  structural  elements,  e.g., 
beams,  the  continuum  (3D  solid)  elements  do  not  include  a  built 
in  kinematic  assumption.  This  characteristic  makes  them  suitable 
for  consistent  incorporation  of  other  effects,  such  as  thermomechanical 
coupling  in  future  studies.  Note  that  in  the  current  study,  no 
coupling,  and  no  material  or  boundary  nonlinearities  are  included. 

Conventional  approaches  in  the  stability  analysis  of  structures 
often  make  use  of  static  considerations,  but  in  fact,  even  when 
the  loads  are  applied  statically,  buckling  and  snap-through  are 
inherently  dynamic  processes  and  a  full  description  of  the 
structural  behavior  can  be  obtained  only  through  a  dynamic 
analysis  [3].  Numerical  simulations  of  such  phenomena  require 
access  to  stable  time-stepping  schemes  and  in  general  to  robust 
simulation  environments. 

Unfortunately,  due  to  the  kinematic  assumptions  incorporated 
in  the  structural  elements,  the  use  of  these  elements  coupled  with 
the  time-stepping  integrators  is  prone  to  numerical  difficulties  that 
affect  the  accuracy  of  the  results  as  shown  in  Section  4.  Numerical 
instabilities  mask  the  true  physical  behavior  rendering  the 
structural  response  prediction  inaccurate. 

An  important  component  of  the  simulation  environment  is  the 
time-stepping  scheme.  We  examine  here  the  common  choices  for 
structural  mechanics  simulations:  (1)  the  traditional  Newmark 
integrator  [8]  and  (2)  energy-momentum  conserving  algorithms. 
The  Newmark  method  is  the  most  widely  used  time  integrator  in 
the  area  of  structural  analysis.  For  certain  combinations  of  para¬ 
meters,  it  is  unconditionally  stable  for  linear  problems.  However, 
its  stability  is  not  guaranteed  for  nonlinear  problems.  The  tradi¬ 
tional  time-stepping  algorithms  developed  for  structural  dynamics 
applications  usually  perform  well  for  linear  problems.  In  the  non¬ 
linear  regime,  however,  numerical  instabilities  appear  due  to  the 
energy  increase  of  the  discrete  system.  Hence,  energy-momentum 
schemes  were  developed  to  overcome  the  lack  of  conservation  [9]. 
These  schemes  belong  to  a  class  of  algorithms  designed  to  satisfy 
various  conservation  laws  by  construction.  The  energy-momentum 
algorithms  used  are  (1)  the  algorithm  based  on  the  work  by 
Simo  and  Tarnow  [10],  Simo  et  al.  [11],  and  Gonzalez  [12]  and 
(2)  a  composite  algorithm  based  on  the  trapezoidal  rule  and  the 
three  point  backward  Euler  proposed  by  Bathe  [13]  that  is  stable 
for  large  time  steps.  These  algorithms  will  be  referred  as 
conserving  A  and  conserving  B,  respectively,  in  the  rest  of  this 
paper.  Many  authors  used  conserving  A  to  solve  various  types  of 
problems  [14].  In  a  recent  work,  Garikipati  et  al.  [15]  used  this 
algorithm  to  model  growth  in  biological  tissue.  Bathe  [13]  showed 
that  the  conserving  B  algorithm,  which  is  incorporated  in  the 


commercial  FEM  software  AD1NA  8.7,  is  able  to  solve  a  specific  type 
of  problem  where  the  Newmark  algorithm  is  unstable  and  does  not 
conserve  energy  and  momentum.  Although  the  conserving  B 
scheme  allows  us  to  use  larger  time  steps,  interactions  between 
the  time-stepping  schemes  and  the  spatial  discretizations  might 
still  occur  when  structural  elements  are  used. 

This  paper  is  organized  as  follows.  Section  2  describes  the 
representative  structure  used  throughout  the  paper.  Section  3 
compares  the  results  of  the  static  analyses  obtained  using  different 
element  types  and  formulations.  Section  4  analyzes  the  numerical 
results  of  the  nonlinear  dynamic  simulations  of  a  curved  beam 
undergoing  snap-through.  Various  time  integrators  and  element 
types  and  formulations  are  used.  A  summary  of  our  findings  is 
presented  in  the  concluding  section. 


2.  Representative  structure 

The  representative  system  discussed  in  this  paper  is  a  curved 
beam  that  undergoes  snap-through  under  a  concentrated  load. 
The  geometry  of  the  beam  is  shown  in  Fig.  2.  The  beam  is 
symmetrical  with  an  angle  0  =  5.674°  at  the  supports.  The  thicker 
lines  in  the  figure  represent  the  parts  of  the  beam  that  are 
clamped;  their  horizontal  length  is  Lc.  The  free  portion  of  the  beam 
has  a  horizontal  projection  Lh.  The  total  horizontal  projection  of  the 
beam,  including  the  clamped  parts  is  L  =  362.204  mm.  The  beam 
can  be  considered  as  a  thin  beam  and  has  a  rectangular  cross 
section  with  depth  d  much  larger  than  the  thickness  t.  Fig.  3  shows 
the  3D  representation  of  the  beam. 

The  dimensions  of  the  beam  are  listed  in  Table  1  and  Table  2 
lists  the  material  properties  of  the  structure.  Without  loss  of 
generality,  the  type  of  load  used  throughout  this  paper  is  a  point 
load  applied  at  the  center  of  the  beam. 

3.  Static  analysis 

In  this  section  we  study  the  effect  of  the  type  of  elements  and 
formulations  on  the  accuracy  of  the  results  in  quasistatic 
simulations.  Based  on  these  preliminary  results,  the  options  that 
are  less  accurate  in  the  static  analysis  are  eliminated,  and  as  a 
result  a  decision  regarding  the  best  elements  to  be  used  in  the 
dynamic  analysis  can  be  made. 

All  simulations  presented  later  in  this  paper  are  performed  with 
meshes  that  ensure  spatial  convergence.  Fig.  4  shows  a  mesh 
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Fig.  2.  Geometry  of  the  structure.  A  curved  beam  clamped  at  the  supports. 
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Fig.  3.  3D  representation  of  the  geometry  of  the  structure  (undeformed  and  post¬ 
snap  configurations). 


Table  1 

Normalized  dimensions  of  the  curved  beam. 


Normalized  dimensions 

Values 

Bottom  radius /L 

5.058 

Free  horizontal  length  (Lh)/L 

0.8415 

Clamp  horizontal  length  (Lc)/L 

0.07925 

Depth  (d)/L 

0.035 

Thickness  ( t)/L 

0.0014 

Rise-to-thickness  (H/t) 

17.67 

Table  2 

Material  properties. 

Properties 

Values 

Young’s  modulus  (N/mm2) 

206,843 

Poisson’s  ratio 

0.28 

Density  (N  s2/mm4) 

7.83  x  10-9 

Number  of  elements  along  the  arch 
Fig.  4.  Mesh  refinement  study  for  3D  solid  quadratic  elements. 


convergence  study  for  3D  solid  quadratic  elements  (27  nodes).  The 
reference  number  of  elements  is  taken  to  be  3072  elements  along 
the  arch.  The  coarsest  mesh  utilized  has  6  elements  along  the  arch 
and  gives  a  relative  error  of  98.2%  in  estimating  the  displacement 
after  snap-through.  Based  on  this  study,  we  use  a  mesh  of  200 
elements  along  the  arch  to  ensure  less  than  1%  relative  error.  A 
similar  study  was  performed  for  discretization  with  beam 
elements  [16]. 

The  structure  is  symmetrical,  but  since  nonsymmetrical 
configurations  after  bifurcation  exist  under  certain  conditions,  all 
simulations  were  performed  on  the  full  beam. 


3.1.  Incremental  loading  using  the  Newton-Raphson  scheme 

The  first  static  analysis  is  performed  with  a  load  control 
algorithm  and  uses  the  Newton-Raphson  iterative  method.  An 
increasing  point  load  is  applied  at  the  midspan  of  the  beam. 

The  load-deflection  curves  for  the  static  analyses  show  that  the 
2D  beam  elements  with  and  without  shear  deformation,  the  3D 
solid  linear  elements  with  enhanced  formulation,  and  the  3D  solid 
quadratic  elements  estimate  similar  snap-through  loads  (Fig.  5), 
Pent  =  1  .51  N.  However,  the  3D  solid  linear  elements  with  displace¬ 
ment  and  mixed  (B-bar)  formulations  experience  locking; 
enhanced  strain  formulations  can  be  used  instead  to  avoid  this 
problem  while  keeping  the  elements  linear,  but  the  use  of  higher 
order  elements  provides  the  same  algorithmic  improvement 
without  the  added  computational  burden  (the  enhanced  elements 
require  a  static  condensation  step  and  inversion  in  every  element) 
and  without  the  increased  storage  requirements  (additional  local 
variables  inside  every  element).  Therefore,  the  2D  beam  elements 
and  3D  solid  quadratic  (27  nodes)  elements  will  be  the  elements 
used  in  the  transient  analysis.  Note  that  there  is  no  warping 
observed  in  the  simulation  that  uses  3D  solid  quadratic  elements. 
Therefore  the  use  of  beam  elements  is  acceptable  even  though  they 
cannot  capture  warping.  Also  note  that  no  torsional  deformation  is 
observed  when  the  discretization  uses  3D  solid  quadratic 
elements. 

We  compared  the  values  of  the  buckling  load  and  the  displace¬ 
ment  after  buckling  obtained  from  the  static  analysis  with  the 
analytical  approximation  proposed  in  Bradford  et  al.  [4],  They  have 
shown  that  these  values  provide  reasonable  approximations  for 
the  symmetric  buckling  of  fixed  arches.  From  the  comparison,  we 
obtained  a  relative  difference  of  4.93%  for  the  buckling  load  and 
a  relative  difference  of  5.47%  for  the  displacement  after  buckling. 

Simpler  models  (link  systems)  were  also  used  for  qualitative 
comparisons  [17]  and  an  extensive  investigation  of  the  post 
snap-through  behavior  of  the  curved  beam  including  identification 
of  solutions  with  non-symmetric  deformation  response  and 
transient  responses  under  a  variety  of  forcing  patterns  is  presented 
in  [18]. 

3.2.  Hysteresis  analysis 

The  hysteresis  analysis  illustrates  the  path-dependence 
behavior.  The  structure  is  loaded  with  a  concentrated  force  (at 
the  middle  of  the  beam)  increased  up  to  a  value  slightly  larger  than 
the  snap-through  value  and  then  it  is  unloaded  with  the  same 
rate.  The  results  under  a  load  controlled  scenario  confirm  that  the 


Fig.  5.  Load-displacement  curve  from  static  analysis  with  Newton-Raphson. 
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structure  experiences  hysteresis  causing  the  equilibrium  path  to  be 
dependent  on  the  history  (Fig.  6). 

3.3.  Path  following  analysis  using  an  arclength  method 

In  order  to  obtain  the  whole  equilibrium  path,  including  the 
unstable  configurations,  an  arclength  algorithm  based  on 
Schweizerhof  and  Wriggers  [19]  is  used.  The  method  controls 
neither  the  load  nor  the  displacement  but  traces  the  equilibrium 
configurations,  both  stable  and  unstable  (Fig.  7),  with  increments 
in  the  length  of  an  arc  on  this  curve.  Even  though  unstable  paths 
are  hard  to  capture  experimentally,  having  access  to  the  whole 
solution  space  provides  significant  insight  into  the  system 
behavior.  If  incremental  loading  is  used,  the  system  advances  to 
point  1  then  snaps  through  to  point  2.  If  the  load  is  increased,  it 
will  travel  to  point  3.  When  the  unloading  starts,  the  system  will 
travel  back  from  point  3  to  point  2  then  to  point  4  before  snapping 
back  to  point  5  and  then  returning  to  the  original  (undeformed) 
position. 

Static  analyses  provide  some  information,  such  as  the  load  level 
at  which  snap-through  and  snap-back  occur  for  specific  material 
properties  and  geometry.  However,  a  dynamic  analysis  is  required 
to  capture  the  transient  behavior  after  the  snap-through  occurs. 
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Fig.  6.  Load-displacement  curve  from  hysteresis  analysis. 


Fig.  7.  Load-displacement  curve  from  static  analysis  with  arclength  method. 
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Fig.  8.  Load  for  transient  analysis. 

4.  Transient  analysis 

This  section  analyzes  the  performance  of  time  integrators  when 
applied  to  transient  problems  involving  snap-through,  i.e.,  the 
dynamic  jump  of  the  system  from  a  “quasi  static”  configuration 
to  oscillations  about  a  remote  equilibrium  configuration. 

The  concentrated  force  applied  at  the  midspan  of  the  curved 
beam  ramps-up  to  some  value  above  the  snap-through  load  and 
then  remains  constant  (Fig.  8). 

In  the  following  simulations,  the  choice  of  finite  elements  used 
to  discretize  the  system  is  based  on  the  preliminary  results  of  the 
previous  section:  2D  beam  elements  with  and  without  shear 
deformation  and  3D  solid  quadratic  (27  nodes)  elements.1 

Due  to  the  snap-through  phenomenon,  the  structure  undergoes 
changes  in  the  stability  behavior.  The  study  of  the  accuracy  and 
robustness  of  the  numerical  methods  for  transient  simulations  of 
structures  that  are  likely  to  traverse  unstable  regimes  raises  the 
following  issues:  (1)  the  existence  of  a  critical  time  step  for  the 
stability  of  the  numerical  integrator,  (2)  the  introduction  of 
unwanted  artificial  numerical  damping,  and  (3)  the  lack  of  exact 
energy  conservation  and  other  numerical  pathologies.  These  issues 
are  discussed  in  detail  in  the  rest  of  this  section. 

To  quantify  the  robustness  of  the  numerical  method,  one  of  the 
measures  used  is  the  total  energy.  The  values  of  the  algorithmic 
parameters  utilized  for  the  Newmark  and  the  energy-momentum 
schemes  ensure  energy  conservation  for  linear  problems.  We  will 
use  the  deviation  from  this  norm  as  a  measure  of  the  loss  of 
numerical  robustness  and  of  the  increase  in  the  likelihood  of  the 
numerical  nature  of  the  instability  (if  instability  is  encountered). 

4.1.  Critical  time  step 

We  first  present  the  numerical  studies  performed  in  order  to 
identify  the  critical  time  steps  for  both  the  Newmark  method 
and  the  conserving  A  algorithm  for  different  finite  element 
formulations  (2D  beams  and  3D  solids).  The  critical  time  steps 
are  identified  by  systematic  numerical  experimentations.  We  also 
present  simulation  results  obtained  using  the  conserving  B 
algorithm,  where  we  will  show  that  there  is  no  critical  time  step 
(the  simulation  does  not  “blow  up”)  but  accurate  results  can  be 
obtained  only  for  small  time  steps. 

Simulations  with  2D  beam  elements  without  shear  may  con¬ 
verge  even  for  At  larger  than  the  critical  time  step.  However, 
for  each  At,  a  different  structural  behavior  is  identified.  This 


1  All  elements  account  for  large  deformation. 
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observation  clearly  suggests  that  even  though  the  nonlinear  solver 
converges,  this  is  not  necessarily  a  solution  of  the  physical  system  but 
rather  a  numerical  artifact.  An  example  of  such  unstable  solution 
(At=10~2s)  suggests  that  the  structure  undergoes  intermittent 
snap-through  before  it  settles  into  periodic  oscillations  around  an 
intermediate  point  that  is  neither  the  original  configuration  nor  the 
snap-through  configuration  (Fig.  9(a)).  The  total  energy  plot  for  this 
simulation  shows  that  the  energy  is  not  conserved;  it  increases 
significantly  when  the  response  settles  into  oscillations  about  this 
intermediate  configuration  (Fig.  9(b)). 

The  simulations  using  the  2D  beam  elements  without  shear 
deformations  are  repeated  with  the  conserving  A  algorithm 
(Fig.  10(a)).  The  energy  plot  for  an  unstable  time  step  also  shows 
that  the  energy  increases  rapidly  at  a  certain  point  in  the  simula¬ 
tion  that  corresponds  to  the  time  increment  where  the  response 
becomes  unphysical  (Fig.  10(b)).  The  total  energy  obtained  in 
simulations  that  converge  to  nonphysical  solutions  is  much  larger 
when  the  conserving  A  scheme  is  used  than  when  using  the 
Newmark  method.  This  makes  the  identification  of  the  unphysical 
solutions  much  easier. 

The  result  of  a  stable  simulation  shows  that  the  structure  expe¬ 
riences  periodic  oscillations,  which  begin  right  after  snap-through 
occurs  (Fig.  11(a)).  Even  though  a  stable  numerical  behavior  is 
exhibited,  we  observe  that  the  numerical  analysis  provides  a 
damped  response.  This  response  is  obtained  from  the  Newmark 
method  with  a  combination  of  parameters  that  should  have 
precluded  numerical  dissipation,  so  clearly  this  damping  is  a 
numerical  artifact,  as  will  be  discussed  in  more  detail  in  the  next 
subsection. 
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Fig.  9.  2D  beam  without  shear  deformation;  Newmark  method;  At  =10  2s 
(a)  displacement,  (b)  total  energy. 
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Fig.  10.  2D  beam  without  shear  deformation;  conserving  A  scheme;  At  =10  2s 
(a)  displacement,  (b)  total  energy. 


The  transient  simulation  performed  with  3D  solid  quadratic 
elements  and  a  stable  time  step  shows  periodic  oscillations  in 
the  post-snap  response  of  the  structure;  the  amplitude  of  oscilla¬ 
tions  stays  constant  in  the  region  of  constant  load  (Fig.  11(b)), 
not  showing  the  dissipation  that  appears  when  the  2D  beam 
elements  without  shear  deformation  are  used  in  the  simulation. 

The  critical  time  step  obtained  from  systematic  numerical 
experimentations  is  summarized  in  Table  3.  The  preceding  study 
shows  that  the  critical  time  step  depends  on  the  type  of  elements 
used  to  discretize  the  structure.  Studies  performed  using  other 
loading  patterns  show  that  these  bounds  also  depend  on  the  load 
pattern  used  (results  are  not  shown  here). 

Simulations  performed  using  the  conserving  B  scheme  [20]  are 
numerically  stable  at  large  time  steps.  However,  this  scheme  intro¬ 
duces  numerical  damping,  which  affects  the  accuracy  of  the  solu¬ 
tions  greatly  when  a  large  time  step  is  used  (Fig.  12(a)).  When 
small  time  steps  are  used  (Fig.  12(b)),  the  solution  matches  the 
results  obtained  using  the  Newmark  method  and  the  conserving 
A  algorithm. 

4.2.  Unwanted  artificial  numerical  damping 

The  discussion  presented  in  the  previous  section  shows  that  the 
2D  beam  element  without  shear  deformation  introduces  artificial 
numerical  dissipation  in  the  transient  analysis  of  structures  under¬ 
going  changes  in  stability  for  time  steps  smaller  than  the  Atcr.  This 
behavior  is  observed  even  though  the  algorithm  used  is  supposed 
to  conserve  energy,  which  suggests  that  there  exists  some 
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Fig.  11.  Displacement  vs.  time;  Newmark  method  (a)  2D  beam 


Time  [s] 

no  shear  (At  =  5  ■  1 0~5  s),  (b)  3D  solid  quadratic  ( At  =  1 0~4  s). 


Table  3 

Bounds  of  the  critical  time  step  for  different  finite  elements. 

Element  type  and  formulation 

Critical  time  step 

2D  beam  element  with  shear  deformation 

2D  beam  element  without  shear  deformation 

3D  solid  quadratic  (27  nodes)  element 

A t„  «  10“4  S 

Atcr«5  10~5  s 

A  ter  «  10~4  S 

correlation  between  the  kinematic  assumptions  built  into  a 
“structural”  finite  element  and  the  numerical  stability  and 
conservation  properties  in  the  transient  simulation  when  the 
problem  is  highly  nonlinear.  (Note  that  the  energy  conservation 
property  is  proven  to  hold  in  the  discretized  system  for  linear 
problems,  but  does  not  hold  for  nonlinear  problems.) 

The  transient  analyses  using  2D  beam  elements  without  shear 
deformation  exhibit  more  numerical  damping  for  smaller  time 
steps  for  both  the  Newmark  method  and  the  conserving  A  algo¬ 
rithm  (Fig.  13).  The  responses  obtained  with  the  Newmark  method 
show  more  dissipation  than  those  obtained  with  the  conserving  A 
scheme.  Therefore,  in  addition  to  the  time  step  size,  the  amount  of 
dissipation  also  depends  on  the  time  integrator  used  in  the  analysis. 
The  relationship  between  dissipation  and  log  At  is  approximately 
linear  (Fig.  13  with  log  scale  representation  on  the  horizontal  axis). 
The  dissipation  might  be  due  to  the  fact  that  the  use  of  no-shear 
beam  formulations  for  dynamic  problems  leads  to  a  discretized 
system  that  is  parabolic  and  consequently  inadequate  for  wave 


propagation  studies  [21  ].  It  is  therefore  expected  that  the 
Newmark  method  that  is  designed  for  hyperbolic  and  hyperbolic- 
parabolic  problems  may  encounter  difficulties  in  the  parabolic 
case. 

To  isolate  the  possible  cause  of  the  artificial  numerical  damping 
in  analyses  with  2D  beam  elements  with  shear  deformation, 
simulation  of  one  other  transient  problem  is  performed;  a  curved 
cantilever  beam,  essentially  representing  half  of  the  representative 
structure  described  in  Fig.  2  that  has  the  same  properties  and  is 
loaded  at  the  free  end  with  a  concentrated  force  increased  in  the 
first  two  seconds  and  then  decreased  to  zero  rapidly  (Fig.  14). 
Notice  that  snap-through  is  not  exhibited  here,  reducing  the 
problem  to  nonlinear  deformation/vibration  only.  In  this  problem, 
the  beam  is  discretized  with  2D  beam  elements  (with  and  without 
shear  deformation). 

When  the  beam  is  discretized  with  2D  beam  elements  without 
shear  deformation,  damping  is  also  present  in  the  case  of  the 
curved  cantilever  beam.  However,  the  damping  is  only  noticeable 
when  the  simulation  time  is  long.  Therefore,  we  can  conclude  that 
snap-through  aggravates  the  amount  of  the  damping. 

Comparison  of  the  responses  for  At=  1CT4  s  obtained  with  2D 
beam  elements  with  and  without  shear  deformation  shows  that 
the  beam  element  with  small  rotation  assumption  and  no  shear 
is  clearly  more  rigid  than  the  beam  element  with  large  rotation 
and  shear,  both  in  the  initial  quasi-static  loading  phase  and  in 
the  amplitudes  of  the  oscillations  in  the  transient  regime 
(Fig.  15).  Since  the  beam  is  very  thin  and  flexible,  the  inclusion  of 
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Fig.  12.  Displacement  vs.  time;  3D  solid  quadratic;  conserving  B  method  (a)  At=  10  3  s,  (b)  At=  10  4s. 
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Fig.  13.  Decrease  in  amplitude  in  the  simulations  using  2D  beam  elements  without 
shear  deformation. 


Fig.  14.  Load  for  transient  analysis  of  the  curved  cantilever  beam. 


the  shear  effects  should  not  bring  significant  changes  to  the  results 
and  we  suspect  that  the  cause  of  this  more  rigid  behavior  is  the 
limiting  assumption  of  small  rotation. 


4.3.  Algorithmic  conservation  of  energy 

In  addition  to  the  existence  of  a  critical  time  step  and  the  intro¬ 
duction  of  artificial  numerical  damping,  the  transient  analysis  of  a 
curved  beam  experiencing  snap-through  also  raises  several  issues 
concerning  the  conservation  of  energy  of  the  system. 

The  analysis  using  2D  beam  elements  shows  that  the  energy 
conservation  is  not  satisfied  throughout  the  simulation.  The  total 
energy  plot  for  2D  beam  elements  without  shear  deformation 
shows  that  the  total  energy  is  significantly  larger  (Fig.  16(a))  than 
the  total  energy  in  the  analysis  with  3D  solid  quadratic  elements 
(Fig.  16(b)).  The  total  energy  includes  the  kinetic  energy,  strain 
energy,  and  the  work  applied.  Under  the  energy  conservation 
measure,  the  conserving  A  algorithm  performs  better  than  the 
Newmark  method.  Simulations  performed  using  the  conserving  B 
algorithm  clearly  do  not  satisfy  the  conservation  of  energy; 
artificial  numerical  damping  is  introduced  into  the  system  [20]: 
the  larger  the  time  step,  the  larger  the  decrease  in  energy  (Fig.  17). 

Pathologies  in  the  variation  of  energy  also  appear  when  2D 
beam  elements  with  shear  deformation  are  used  for  the  spatial 
discretization  of  the  system.  The  results  obtained  from  analyses 
performed  with  the  Newmark  method  (At=  10~4s)  for  a  system 
discretized  with  2D  beam  elements  with  shear  deformation  show 
that  the  refinement  of  the  mesh  increases  the  energy  of  the  system 
(Fig.  18).  Note  that  the  energy  does  not  change  when  a  smaller 
time  step  is  used,  showing  that  the  increase  in  energy  depends 
only  on  the  mesh  refinement. 

Further  examination  shows  that  as  the  mesh  is  refined,  it  is  the 
kinetic  energy  that  changes  significantly  in  magnitude  while  the 
strain  energy  remains  the  same.  This  phenomenon  is  observed  only 
when  the  2D  beam  elements  with  shear  deformation  capabilities 
are  used  in  the  analysis.  The  analyses  using  2D  beam  elements 
without  shear  deformation  and  3D  solid  quadratic  elements  show 
no  dependence  of  the  kinetic  energy  on  the  level  of  mesh 
refinement. 

Several  attempts  were  made  to  understand  this  behavior.  A 
closed-form  solution  based  on  the  assumptions  of  negligible 
displacement  in  the  X  direction  and  small  rotation  shows  that 
the  kinetic  energy  for  one  element  with  length  h  and  two  elements 
with  length  /t/2  in  one  time  step  is  the  same,  but  once  again,  such 
results  no  longer  hold  for  large  deformation  [16]. 

Simulations  performed  on  the  curved  cantilever  beam  used  in 
the  previous  subsection  examined  whether  the  increase  in  energy 
appears  in  other  transient  problems  solved  with  2D  beam  elements 
with  shear  deformation  capabilities.  Recall  that  this  formulation 
also  accounts  for  large  deformation  and  large  rotation.  The  geo¬ 
metry  represents  half  of  the  representative  structure  described  in 
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Fig.  15.  Curved  cantilever  beam;  Newmark;  At=  10  4s;  various  2D  beam  elements  (a)  displacement  vs.  time,  (b)  zoom  in. 
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Fig.  17.  Total  energy  vs.  time;  3D  solid  quadratic;  conserving  B  method  (a)  At=  10  3  s,  (b)  At=  10  4  s. 
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Fig.  18.  Arch  fixed  at  both  ends;  2D  beam  with  shear;  Newmark 
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;  At=  10  4  s;  various  meshes  (a)  total  energy  vs.  time,  (b)  zoom  in. 


Fig.  2  and  the  loading  is  shown  in  Fig.  14.  The  kinetic  energy  plots 
show  that  there  is  an  increase  in  the  kinetic  energy  when  the  mesh 
is  refined  (Fig.  19)  for  this  case  as  well.  Therefore,  we  can  conclude 
that  the  dependency  of  the  kinetic  energy  on  the  mesh  refinement 
is  a  feature  of  the  2D  beam  with  shear  elements  and  not  necessar¬ 
ily  induced  by  the  type  of  the  problem  (e.g.,  snap-through). 


5.  Conclusions 

This  paper  analyzes  the  performance  of  several  finite  element 
formulations  (2D  and  3D)  and  the  stability  of  the  time-stepping 
schemes  in  simulating  a  curved  beam  undergoing  snap-through 
by  identifying  the  important  features  that  affect  the  numerical 
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Fig.  19.  Curved  cantilever  beam;  Newmark;  2D  with  shear;  A t=  10  4s;  various  meshes  (a)  kinetic  energy  vs.  time,  (b)  zoom  in. 


accuracy  and  robustness  and  the  region  where  the  schemes  are 
stable  for  such  simulations.  We  examine  the  interaction  between 
the  two  most  important  components  of  the  finite  element  analysis 
applied  to  structural  dynamics  problems;  (1)  the  time-stepping 
schemes  and  (2)  the  finite  element  formulations  used  to  spatially 
discretize  the  structure.  The  integrators  studied  are  (1 )  the  Newmark 
method,  (2)  a  conserving  integrator  (referred  as  conserving  A) 
[10-12],  and  (3)  a  composite  integrator  (referred  as  conserving 
B)  proposed  by  Bathe  [13],  The  finite  element  formulations  used 
are  the  2D  beam  elements  with  and  without  shear  deformation 
capabilities  and  the  3D  solid  elements. 

The  study  shows  that  the  Newmark  method  and  the  conserving 
A  scheme  have  a  restrictive  bound  on  the  size  of  the  time  step  to 
ensure  numerical  stability.  The  critical  time  step  depends  on  the 
finite  element  formulations  used  to  discretize  the  structure  and 
on  the  loading  pattern.  Simulations  performed  using  the 
conserving  B  algorithm  do  not  have  this  restriction.  However  the 
numerical  damping  introduced  by  this  algorithm  increases  with 
the  increase  in  the  time  step  and  greatly  affects  the  accuracy  of 
the  solution  when  large  time  steps  are  used. 

The  study  also  shows  that  the  spatial  and  temporal  discretiza¬ 
tions  may  interact  and  such  interactions  may  induce  unwanted 
numerical  effects  such  as  artificial  damping,  lack  of  energy  conser¬ 
vation,  and  most  importantly,  misleading  numerical  results  that 
seem  to  indicate  a  chaotic  response  when  in  fact  the  simulation 
simply  converged  to  non-physical  solutions.  These  issues  are  very 
severe  when  structural  elements  are  used  to  discretize  the  beam 
while  the  use  of  3D  solid  elements  has  less  severe  effects.  There¬ 
fore,  we  recommend  the  use  of  3D  solid  elements  in  solving 
snap-through  problems. 

Due  to  the  nonlinearity  of  the  problem,  the  conservation  of 
energy  is  sometimes  lost  when  structural  elements  (e.g.,  beams) 
introduce  kinematic  assumptions  in  the  discretization  of  the  struc¬ 
ture.  The  total  energy  is  approximately  conserved  when  3D  solid 
quadratic  elements  are  used  with  the  Newmark  method  and 
conserving  A  algorithm.  When  the  conserving  B  algorithm  is  used, 
the  conservation  of  energy  is  not  satisfied  because  numerical 
damping  is  introduced  to  the  system.  The  inspection  of  energy 
plots  for  the  analyses  performed  with  beam  formulation  with  shear 
deformation  also  shows  another  numerical  artifact:  the  kinetic 
energy  is  mesh  dependent  (increases  as  the  number  of  elements 
increases;  reducing  the  time  step  does  not  eliminate  this  unwanted 
behavior). 

In  conclusion,  we  have  shown  that  several  of  the  currently 
available  finite  element  formulations  are  not  robust  or  accurate 
enough  to  simulate  snap-through.  In  particular,  assumptions 


built  into  the  formulation  of  structural  elements  (e.g.,  beam 
elements)  lead  to  unwanted  numerical  behavior  that  is  often  times 
amplified  by  the  presence  of  the  snap-through  phenomenon  in  the 
system. 
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